{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Earth orbit\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "year"
      ],
      "text/latex": [
       "$\\mathrm{year}$"
      ],
      "text/plain": [
       "<Unit('year')>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Here are the units we'll need\n",
    "\n",
    "s = UNITS.second\n",
    "N = UNITS.newton\n",
    "kg = UNITS.kilogram\n",
    "m = UNITS.meter\n",
    "year = UNITS.year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>R</th>\n",
       "      <td>[147000000000.0 meter, 0.0 meter]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>V</th>\n",
       "      <td>[0.0 meter / second, -30300.0 meter / second]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "R                [147000000000.0 meter, 0.0 meter]\n",
       "V    [0.0 meter / second, -30300.0 meter / second]\n",
       "dtype: object"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# And an inition condition (with everything in SI units)\n",
    "\n",
    "# distance from the sun to the earth at perihelion\n",
    "R = Vector(147e9, 0) * m\n",
    "\n",
    "# initial velocity\n",
    "V = Vector(0, -30300) * m/s\n",
    "\n",
    "init = State(R=R, V=V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>R                [147000000000.0 meter, 0.0 me...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>G</th>\n",
       "      <td>6.674e-11 meter ** 2 * newton / kilogram ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>m1</th>\n",
       "      <td>1.989e+30 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>r_final</th>\n",
       "      <td>701879000.0 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>m2</th>\n",
       "      <td>5.972e+24 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>31556925.9747 second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "init       R                [147000000000.0 meter, 0.0 me...\n",
       "G              6.674e-11 meter ** 2 * newton / kilogram ** 2\n",
       "m1                                        1.989e+30 kilogram\n",
       "r_final                                    701879000.0 meter\n",
       "m2                                        5.972e+24 kilogram\n",
       "t_end                                   31556925.9747 second\n",
       "dtype: object"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Making a system object\n",
    "\n",
    "r_earth = 6.371e6 * m\n",
    "r_sun = 695.508e6 * m\n",
    "\n",
    "t_end = (1 * year).to_base_units()\n",
    "\n",
    "system = System(init=init,\n",
    "                G=6.674e-11 * N / kg**2 * m**2,\n",
    "                m1=1.989e30 * kg,\n",
    "                r_final=r_sun + r_earth,\n",
    "                m2=5.972e24 * kg,\n",
    "                t_end=t_end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here's a function that computes the force of gravity\n",
    "\n",
    "def universal_gravitation(state, system):\n",
    "    \"\"\"Computes gravitational force.\n",
    "    \n",
    "    state: State object with distance r\n",
    "    system: System object with m1, m2, and G\n",
    "    \n",
    "    returns: Vector\n",
    "    \"\"\"\n",
    "    R, V = state\n",
    "    G, m1, m2 = system.G, system.m1, system.m2\n",
    "        \n",
    "    # make sure the result is a vector, either\n",
    "    # by forcing it, or by putting v.hat() on\n",
    "    # the left\n",
    "    \n",
    "    force = -G * m1 * m2 / R.mag2 * R.hat()\n",
    "    return Vector(force)\n",
    "\n",
    "    force = -R.hat() * G * m1 * m2 / R.mag2\n",
    "    return force"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\\[\\begin{pmatrix}-3.6686485997501037×10<sup>22</sup> & -0.0\\end{pmatrix} newton\\]"
      ],
      "text/latex": [
       "$\\begin{pmatrix}-3.6686485997501037\\times 10^{22} & -0.0\\end{pmatrix}\\ \\mathrm{newton}$"
      ],
      "text/plain": [
       "array([-3.6686486e+22, -0.0000000e+00]) <Unit('newton')>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "universal_gravitation(init, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The slope function\n",
    "\n",
    "def slope_func(state, t, system):\n",
    "    \"\"\"Compute derivatives of the state.\n",
    "    \n",
    "    state: position, velocity\n",
    "    t: time\n",
    "    system: System object containing `g`\n",
    "    \n",
    "    returns: derivatives of y and v\n",
    "    \"\"\"\n",
    "    R, V = state\n",
    "\n",
    "    F = universal_gravitation(state, system)\n",
    "    A = F / system.m2\n",
    "    \n",
    "    return V, A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([     0., -30300.]) <Unit('meter / second')>,\n",
       " array([-0.00614308, -0.        ]) <Unit('newton / kilogram')>)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Always test the slope function!\n",
    "\n",
    "slope_func(init, 0, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here's an event function that stops the simulation\n",
    "# before the collision\n",
    "\n",
    "def event_func(state, t, system):\n",
    "    R, V = state\n",
    "    return R.mag - system.r_final"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "146298121000.0 meter"
      ],
      "text/latex": [
       "$146298121000.0\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "146298121000.0 <Unit('meter')>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Always test the event function!\n",
    "\n",
    "event_func(init, 0, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>Success</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "message    Success\n",
       "dtype: object"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Finally we can run the simulation\n",
    "\n",
    "system.set(dt=system.t_end/1000)\n",
    "results, details = run_euler(system, slope_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_trajectory(R):\n",
    "    x = R.extract('x') / 1e9\n",
    "    y = R.extract('y') / 1e9\n",
    "\n",
    "    plot(x, y)\n",
    "    plot(0, 0, 'yo')\n",
    "\n",
    "    decorate(xlabel='x distance (million km)',\n",
    "             ylabel='y distance (million km)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3xb5bnA8Z8k773tOLZjZ/hNnL13AoQV9iibsgu0XO4tLbTcXspqKS0dlEJbSqGMskdZAcLIDtnDiRPHb2wnjvfe25Z0/ziycYLjKMGyZPv5fj7nI/kcyecJJHr0ruc12e12hBBCCE9jdncAQgghRG8kQQkhhPBIkqCEEEJ4JElQQgghPJKXuwNwN6WULzAbKAGsbg5HCCGGGwswAtiutW7reWHYJyiM5LTB3UEIIcQwtxjY2POEJCij5cRrr71GXFycu2MRQohhpbS0lOuuuw4cn8U9SYJydOvFxcWRkJDg7liEEGK4+tYQi0ySEEII4ZEkQQkhhPBIkqCEEEJ4JElQQgghPJIkKCGEEB5JEpQQQgiPJNPMhfBwVquNlnYrbe2dtLR10tFpw2qzY7PZsVrt2Ox2rDYbVqsdq82O2WTCYjZhdhw9n3tbzPj6WPDz8cLXx4KvtwWz2eTuP6IQvZIEJcQAae+wUl3fSlVdK3WNbTQ0t1Pf1E5DcwcNTe3dPze2tNPS1pWQrHRabS6Ny8fbgp+PcQT5+xAU4E1wgA/BgT4Edz0P8CYk0JeIED8iQv0IDfLFIolNuJgkKCH6gc1mp6quldLqJkormyirbqaqrpXq+lZHUmqhobnjlH632QS+Pl74+1rw9fHC28uMl9mM2WK0jozD3N1SstkdrStHK8tms2O127FZ7XRYrbS1W2l1HO0d3xz1TVBe0+J0TGHBRrKKDPEjIsSP6HB/4iIDiYsMIDYikOAAb0wmSWLi1EmCEsJJdrud6vpWCsoayC9roLiiidIq4yirbjlhS8diNhERanyYhwX5EhLo06Ol4kNIYFdrxQd/X6MLzt/XSEiu+qC32ey0dVhpbe+ktc1KU0sHDc3tjqODRsdjQ3M7tY1t1DgSbl1je3fyzTnO7w7w8yIuIpDYyABGRAaSGBtMYmwQibHBBPh5u+TPI4YWSVBC9KKusY3cojrySxuMhFRaT0FZA02tncd9T1iwL3ERAcRFGh/KUaH+3S2MyFB/QgJ9PG68x2w24e/rhb+vFwQ7/76OThs1DY4WYp3RbVle0+xI2M2UVTfR3NrJoeI6DhXXfev9kaF+JMYGkxQbTGJsMKNHhjJqRAi+3pZ+/NOJwc7jEpRSag6wQmsd4/jZF2gA2nu8bJPW+mzH9SuB32CUa18H3KS1Lh/YqMVgVl3fSm5hLTmFdeQW1pJbVEdlbe9dXUH+3iTFGR+qCTFBji6tQOIiAvDz9bh/Ti7j7WUmJjyAmPCAXq/b7Xbqm9opqzaSVnFlEwWlRsuzqKKRKkdSSz9Y0f0es9lEYkwQo0eGMiYhjNEjQxkdH0qgv7S2hiuP+RellDIBtwJ/OObSZKBaa/2tUuNKqTTgBWA5sAP4HfAmcIZroxWDVUenldzCOg7kVXMgr5qsvGpqGtq+9To/Hwsp8aEkx4d0f8tPig0mLNhXxlWcYDKZCA3yJTTIl9Sk8KOuWW12yqqbKCxrJN/ROj1UVEdBWQNHSo1jzc7C7tcnxAQxITkCNSqCCcnhJMQEe1xLVLiGxyQo4BHgfODXwAM9zs8E0o/znuuBj7XWGwGUUv8L1Cilxmmts10ZrBgcGpvb2XeoigOHjYSUU1hLR+fRY0UBfl6MGRnGmATjm/uYkaHERwfJLDUXsZhNxEcFER8VxJyJ33zvbOuwcqSkvrsVm1tUx5GSegrLGyksb+TLbfkABPp7o0aFMyE5gokpkYxPDsfbS7oGhyJPSlDPaq0fVEqddsz5GUCMUmovEAusB36stS4C0jBaTgBorZuVUgUYrS5JUMNQa3snmYer2ZtdwZ7sCnKL6rDbj35NYmwwE5IjmJAcwfjkcOKjguQbuQfw9baQmhR+VIuro9PG4eKjW7xVda3syipnV5bRk+/jZSYtJZIp46KYOi6aMSNDsVikBsFQ4DEJSmtdfJxLTcDXwKNAB/AX4H1gDhAENB/z+mag945xMeTY7XZyi+rYeaCMPdmVHMirPmo2nZfFhBoVwcTRkUZCGhVOUICPGyMWJ8Pby9ydtC5eMgaAipoWsvKqycyrYl9uFXkl9aRnV5CeXQEcIMDPi0mjo5iWGs3stFjiIgPd+4cQp8xjEtTxaK1/0vNnpdRPgAqlVCJG8vI/5i0BQOMAhSfcoLWtkz3ZFWw/UMb2zDKq61u7r5lMMDYhlCljo5k6Lpq0lIhhNXlhOIgO9yc6fCSLp48EoLahjYycSvbkVLA3p5KSyia2ZZayLbOU5z7IIDE2iFkT4pg9IZYJKRF4Setq0PD4f7lKqUeBN7TWBxynur7+tgKZgOrx2gAgyXFeDCHV9a1szihhe2Ype3MqjxpHigr1Y1ZaHNNTo5k8NopgaSENK2HBviye/k3CKq9uZm9OBTuyytmtyykoa6SgLIf31+YQ6OfFdBXD3EkjmJMWK+uxPJzHJyhgCjBLKXWt4+engE+01hVKqdeBjY5xq83A48BurfVB94Qq+lNVXQub9pbw9d5iMg9XdY8lmUygRoUzOy2WOWlxJI8IkZl1oltMRABnzhnFmXNG0Wm1kXm4iu2ZZew4UEZheSMb9xSzcU8x3l5mZqgYFkyJZ+7EOJnO7oEGQ4K6FWPcKQcj3k+A2wG01hlKqVuAZ4GRwFbgCjfFKfpBdX0rX+8p/lZS6vowmTdpBLMmxBIW7OveQMWg4GUxM2VsNFPGRnPrRZMorWpi2/5SNmWUkHm4iq37S9m6vxQvi4lpqTEsnBLP/MkjJFl5CJP92ClOw4xSKhk4vGrVKhISEtwdzrDU1mFlS0YJq3cUkH6wHFuPpDRzfAwLp46U7hjR77q6jTftLWZfbmX33zsfLzPzJo3g9FmJTE+NlhmBLlZYWMiyZcsAUrTWeT2vDYYWlBiC7HY7B/KqWb2jgA3pRTQ7Sgh5WUzMHh/LommSlIRrRYT4cf7CFM5fmEJNQytb9pWyYXcRGbmVrE8vYn16EWHBvpw2I4EzZiWSEh/q7pCHHUlQYkDVNbbx1bZ8Pt96hJLKpu7zYxPDWDYrkSXTEwgJlEkOYmCFB/uxfH4yy+cnU17dzNpdhazekU9RRRMfrMvlg3W5pMSHcO78ZE6bkSBfnAaIJCjhcl2tpc825bFxT3H3OqWIED9On2l8O02KC3FzlEIYYiICuPLMVK5YNo6D+TWs3lHA+t1FHC6u5+/v7eWlFftZOiOR5fOTGT1SWlWuJAlKuExzawdrdxXy2aY88krqAWMG3qwJsZy3IJkZ42OlnJDwWCaTschbjYrgtosnsTmjhE835bH/UBUrN+excnMeKimc5QuSWTxtJD5Sib3fSYIS/a6qroWPNxxi5ea87u0pwoJ8OWtuEufMSyY2Qgp9iMHF28vCkukJLJmeQH5pPZ9tzmPNjgJ0fg06v4aXVmRy3sIUzluQTGiQzDDtL5KgRL/JK6nn/bU5rN9dSKfVmBI1ITmCCxalMH9yPN5eMhtKDH5JcSHccekUbjwvjQ3pRaz4+jCHiup4/fMs3l11kDNmJ3HxktEkxJzEBluiV5KgxHdit9vJyK3kvdU57HJsw2U2wcKp8Vy6dAxqVISbIxTCNfx8vThr7ijOnJNERm4l76/NZceBMlZuzuPzLXnMSYvj8tPHMSFF/g2cKklQ4pTY7Xb25lTyxhea/YeqAPD1sXDWnCQuXjJGCnSKYcNkMnUvBi4oa+CDdbms2VnQvQh46rgorj5LMWlMlLtDHXQkQYmTYrfb2ZtdyetfZJF5uBqA4ABvLl4yhvMWpkgdPDGsJcYGc/eV07h++Xg+2XiYjzceYk92JXuyK5k8JoprzlZMHiuJylmSoITTMnIqeXXlgaMS0yVLx3LBohRZFyJED+HBfly/fAKXLB3DRxsO8dH6XDJyK8n4eyUTR0dy3TnjJVE5QRKUOKG8knpeWrGfnY4N4oIDvLn0tLGcv1ASkxB9CQrw4dpzxnPxkjF8vPEQH67LZf+hKn7x96+ZOT6GG89PkwoVfZAEJY6roqaF1z4/wOodBdjt4O/rxeWnj+XCxaMlMQlxEgL9vbn6LMVFi0fz8YZDvLcmh51Z5ezS5Zw+M5HrzhlPjCy/+BZJUOJbWto6eWfVQT5cl0t7pw2L2cTyhclcfZaSNR5CfAcBft5cdZbi3PnJvP3VQT7ddLi7UsWFi0dz9Vmp8uWvB0lQopvdbmdDehH/+ng/VXXGLrWLpsbz/fMmEB8V5ObohBg6QoN8+cElk7lw8Whe/SyLdbsLeX9tDmt2FnDT+WmcPjMRs1RZkQQlDIeL63jugwz25RpTxscmhHLHpVMYnyxrOIRwlbjIQO69fiaXLB3DP97fS9aRGv785m4+25TH7ZdOJjUp3N0hupUkqGGupa2T11Zm8fGGXGx2CAn04Ybz0jhrTpJ8gxNigIxNDOOJuxezdlchL63Yj86v4adPrefsuaO4+YI0gobp8g1JUMPYrqxy/vpuOuU1LZhNcMHCFK47d/yw/ccghDuZTCZOn5nI3IlxvP3VQT5cn8sXW4+wLbOU2y+ZzKKp8ZhMw+tLoySoYaiusY0XPtrHmp2FAIyOD+XuK6cxNjHMzZEJIQL8vLnpgoksm53EM++kk3m4mif+vYPVO2L54eVTiAkfPrP9JEENM5v2FvPXd/dQ39SOj5fZWKOxdAxesq21EB4lMTaYx3+0iC+2HuGlFfvZcaCMu55Yzc0XTmT5/ORh0ZqSBDVMNLV08NwHGazeUQDAlLFR3HXFVJmdJ4QHM5tNnDs/mTkT43ju/Qy+3lvM39/by9Z9pfz3VdOIDPV3d4gu5XSCUkqZgfFADGAFSoEcrbXdRbGJfpKRU8mTb+6ioqYFHy8zN184kfMWpMgkCCEGiYgQP+6/cTYb9xTxt3f3sEuXc9fv1/DDy6awZPrIIduaOmGCUkotAf4HOBPoucGJHahRSq0E/qa13uSaEMWp6rTaePWzA/xnbQ52uzFT6CfXzCAxVvapEWIwWjR1JGkpkTz9djo7DpTxh9d2smVfCXddMY0g/6G3wPe4CUopNQ74B5AEvA9cBmQCVYAZiAamAkuAN5VSucAdWuuDrg5anFhlbQtP/HsHB/KqMZtNXHVmKledlSpjTUIMchEhfjx461y+2HqE5z/cx8Y9xWQX1PLzG2YxLnForZvqqwX1KvCo1vqT41wvcBwrlFI/By5xvGfOdwlIKTUHWKG1jnH87AM8A3wPo2vxT1rrx3u8/krgN8AIYB1wk9aOnfOGqR0HyvjT67toaG4nMtSP+66fxcTRke4OSwjRT0wmE+fMS2by2Ch+98oODhXV8bOnN3DzBRO5cPHoIdPl19fX6Xl9JKejaK3tWuv3gbmnGohSyqSUug34Aui5EOcRQAFjgNnAjUqpGxzvSQNeAG4CIoFs4M1TjWGws9rsvPJpJo88v4WG5nZmqBie+slpkpyEGKLio4L4/d2LuWBhCp1WO//8cB+/eWkbjS0d7g6tXxy3BdXb5AelVCDwrWqhWuvq473nJDwCnA/8Gnigx/kbMVpFNRhjXn8A7gBeAa4HPtZab3TE97+O14zTWmd/h1gGnaaWDv7w2k52HCjDbILrzp3A984YJxMhhBjifLwt3HHZFCaNjeLpt3azZV8p+X9exwO3zB30481ODUgopc5SSuUB9UBFj6PS8dgfntVazwR29LhvGEbXXWaP12UBkx3P03pe01o3Y3Q7TmYYKapo5N6/rGfHgTKCA3x49I4FXHlmqiQnIYaRhVPi+fNPTiMlPoTiyiZ++tR6tu0vdXdY34mzI+Z/A7YBZwALehzzHY/fmda6uJfTXYt0mnucawYCelxv5mg9rw95u7LK+emf11FY3siouGD+9OMlTB0X7e6whBBuEBcZyBP/tZiFU+Npaevk1y9u5a0vNXb74FwN5Ow6qBHA+W6YodfkeOy5Gi0AaOxx/diVaj2vD2mfbznC397bg81mZ/7kEdxzzQz8fWXttRDDmZ+vFz///izeHZnNvz87wKsrs8gvbeDH10zH28vi7vBOirMtqI+Bs1wZSG8c406lGJMkuoznm269zJ7XlFIBGNPie3YJDjl2u53XVmbxzDvp2Gx2rjwzlftvmC3JSQgBGLP8rliWygO3zMXf18L69CIefG7zoJs84ewn2r3AXqXUVcAhwNbzotb6lv4OrId/Aw8ppfZidOndCzzluPY6sFEpdRqwGXgc2D2U12J1Wm387d09fLktH7MJ7rx8KsvnJ7s7LCGEB5qTFsfjP1rEoy9sYV9uFT97egMP/2DeoCk4ezJjUN5AAxCIUVGi5+FKDwL7gP3AduA94FkArXUGcIvj50pgInCFi+Nxm/YOK795aRtfbsvHx9vCL26aI8lJCNGnMQlh/P7uJSTGBlFQ1sB9f1nP4eI6d4flFJMzg2dKqWbgzKFYzkgplQwcXrVqFQkJCe4O57ha2zt57F/bSM+uIDjAmwdvm8f4UbLbrRDCOY3N7Tz20jb25VYR5O/No3fM94jKE4WFhSxbtgwgRWud1/Oasy2oQqC1n+MSTmpu7eDhf24hPbuCsCBfHv/RIklOQoiTEhTgw6O3z2fuxDgaWzp44NlNZB6ucndYfXJ2DOonwPNKqV8DucBRI21a6yE9KcGdmlo6ePifm8k6UkNEiB+/vnPBoF98J4RwD28vC/ffOJs/vb6LDY6JE7+8eS5TUz1zaYqzLaiPgGnAu8BujDGhjB6PwgVa2zp55PktZB2pITrcn9/etUiSkxDiO/GymPnpdTNZNjuRtnYrj7ywhT0H+6veQv9yNkGl9HKM7vEo+llHpzEh4kBeNVGhfjz+o0WMiAp0d1hCiCHAYjbx31dOZ/n8ZDo6bfzqxa0e2d3nbBefn9ZaH3tSKeULPAT8ol+jGuasVhu/f3Unuw9WEBrkw6/uXEBsxOCYFiqEGBzMZhN3XjaF9k4rq7YX8MjzW3jszoWMTQxzd2jdnG1BrVNKTep5Qil1NsaC2Dv6PaphzG638/Q76WzOKCHQ35tf3bGAhBjp1hNC9D+z2cTdV05n0dR4mls7efC5TRwpqXd3WN2cTVCvAGuVUjOUUjFKqTeAz4D1GJUdRD9566uDrNpegK+PhYdvm0dKfKi7QxJCDGEWs4mfXDuT2WmxNDR38NA/N1NV1+LusAAnE5TW+mfA74HVwAFgArBYa32z1tozR9cGobU7C3htZRZmE/zs+lmMT5ap5EII1/P2MnP/DbNJS4mgqq6VR57fQnOr+8siOb3/t9b6dxjTzQOBnw3FRbvutP9QFU+9lQ7ArRdPYs7EODdHJIQYTny8LfzfzXOJjwrkcHE9v/v3DqxW24nf6ELHnSShlKoAeiszYQY+Vkp118ro2p5dnJry6mYee3EbnVYbFyxK4aLFY9wdkhBiGAoJ9OGhH8zjvr9sYFdWOc++n8Fd35vqtnj6msV374BFMYy1d1h5/JXtxhbt42O47eJhtdeiEMLDxEcF8ctb5vJ/f/+alZvzGJsQyjnzkt0SS19bvr88kIEMV899kEFOQS0xEQHce91MLLILrhDCzcYnR3DXFdN48o1dPPufDJJHhKDcUF7N6TEo0f++2naEz7ccwdvLzP/eOJvgAB93hySEEACcMSuRCxam0Gm18fjL26lpGPhyrJKg3KSgrIG//8eoEvWjy6cwNsFzFscJIQTALRdNYkKyMbPv9//eidU2sFvHS4Jyg45OG396fSftHVbOmJXImXNGuTskIYT4Fm8vM/ffOJvwYF8yciv5z5rsAb2/JCg3ePNLTU5hHTHh/tx+iUyKEEJ4rogQP3589QwAXluZRXZBzYDd26lafEopf+BOYCbGzrpHjeRrra/s/9CGpgOHq3l31UFMJvjJtTMJ9Pd2d0hCCNGnGeNjuGjxaD7acIg/vLqTp35yGn6+zpZyPXXOtqCeA34N+APNQNMxh3BCR6eVp97ajc0Ol502lomjI90dkhBCOOXG89MYFRdMcWUTL3y8f0Du6WwKPBe4Vmv9oSuDGereXZ1DUUUjI6ODuO5cKWEohBg8fLwt3Hv9LO55ci0rN+exdPpIJo2Jcuk9nW1BdQAHXRnIUFdU0cjbXxn/Ce+6YireXhY3RySEECcneUQI3zsjFYBn3tlDe4fVpfdzNkE9CfxWKeWZ+wJ7OLvdzt/e3UOn1caZs5OY7OJvHUII4SpXnjmOhJgg40v3Kte2W5zt4rsSmAKUKqUagPaeF6UWX9+27Cthb04lIYE+3HzhRHeHI4QQp8zby8J/XTGN+/+6kfdWZ7N0egKJsa7Zs87ZBPWMS+4+DHR02nhxRSYA154znpBAqRYhhBjcJo6O5Jx5o/h8yxH++UEGj9w+H5Op/8u0OZWgetblU0qFAGatdW2/R9MHpdQtwD+Ath6n7wLewEig3wOswJ+01o8PZGx9Wbk5j5LKJkZGB3LOPFmQK4QYGr6/fAIb04vYfbCCXbqcmeNj+/0eTi/UVUr9UClVANQAVUqpEqXU/f0e0fHNAP6otQ7qcbwMPAIoYAwwG7hRKXXDAMZ1XE0tHbzxhQbgpgsm4mWRddFCiKEhNMiXG85PA6CixjU78Dq7UPde4JfAY8BGjIW6C4H7lVItWuunXBLd0WYCvd3nRuAmrXUNUKOU+gNwB8Y29W614utDNDS3M3F0JHNlA0IhxBBz3oIUZk+IIzLUzyW/39kxqLuAO7XWb/Q497VS6gjGAl6XJiillAVjksb3lVJ/wlgs/DxGl98IILPHy7MAt9cPam3r5MN1hwC45mzlkv5ZIYRwt+hwf5f9bmcTVDSwvZfzO4GE/gunz/vvAF4GLgMmAB8CXTMOmnu8thkIGICY+rRySx4Nze2oUeFMGSvTyoUQ4mQ5m6D2AVcAx04+uAqjxeJSWutSYGmPU+lKqaeB5Y6fe6bwAKDR1TH1pb3DyvtrcwC46sxUaT0JIcQpcDZBPQh8opSaD2x2nJuPUQLpMlcE1pNSaiJwpdb6oR6nfYBWoBRjkkSR4/x4ju7yG3Ab9xRTXd9G8ogQZk3o/5ktQggxHDg7zfwLpdQy4G7g+0ALcACYrbXe48L4utQCP1VKFQIvANOB/wb+C9gPPKSU2gsEAffi4jGxE/l002EALlw8WlpPQghxipyul661Xg+sd2Esfd27SCl1EfAERtmlSuBXWut3lVIrgD9iJCozRuX1Z90RJ0BOYS36SA2Bfl4smT7SXWEIIcSgd9wEpZR6G7hNa13veH5cA7EflNZ6NTCrl/OtGLMM73J1DM74bFMeAMvmJOHn4/r9UoQQYqjq6xO0CbD3eC5OoLW9kw3phQAsn5/s3mCEEGKQO26C0lrf3NtzcXzbM8toabOSmhRGQoxriicKIcRw0VcX33lO/g671vqzfopnUFu/22g9LZk+EEvDhBBiaOuri2+Fk7/DDgz73feaWjrYcaAckwkWT5PJEUII8V311cUnlU1Pwu6D5XRabUwcHUlEiGvqUgkhxHDSVxef0+WCtNbNJ37V0LY9swxAFuYKIUQ/6auLr5FvZvEdjwnp4sNms7MrqxyA2ZKghBCiX/SVoM7gxAlKAHkl9dQ2thEV5k9SnMzeE0KI/tDXGNTaAYxjUMs8XAXA5DGRUtpICCH6SV9jUOVAmta6UilVQR+tKa11jCuCGyz2HzISVFpKpJsjEUKIoaOvLr77gAbH83sHIJZBK/NwNQATR0uCEkKI/tJXF9/LvT0XR6tpaKW6vpUAPy8SYoLcHY4QQgwZTlUzVUqFAT8GJgK+x17XWl/Uz3ENGnnF9QAkjwiR8SchhOhHzpbbfh2YCXyKsdWFcMgr+SZBCSGE6D/OJqglwFla680nfOUwc6RUEpQQQriCs+WMDmEsyhXHKKs2imiMiAp0cyRCCDG0ONuC+iHwjFLqaeAwYOt50bHb7rBU7khQMRFOV4YSQgjhBGcT1AxgMvBCL9eGbamjTquNytoWTCaIDvN3dzhCCDGkONvF90vgN0AsEHzMMWwHX+oa27DZITTQF2+vYZmjhRDCZZxtQfkAL2mtK1wZzGDT2NIBQHCgt5sjEUKIocfZFtTfgR8rpWSPqB4am40EFeTv4+ZIhBBi6HG2BTUOuAC4QSmVB3T0vKi1ntO/YQ0Ojc3tAAQFSAtKCCH6m7MJKsNxiB46rMZkRh8ZfxJCiH7nVILSWj/i6kC+C6XUVOBZYArGmq1btNbbXX1fq9Uo8G4xyxIxIYTob8cdU1JKfaWUcrrrTim1UCm1un/Ccp5Sygf4EHgLCAMeA75QSrl8dqHVZiQos0USlPAMVcW7yFj/GDu/uI+M9Y9RVbzL3SEJccr6akHdD/xDKWUF/gOsBDK11u0ASilfYCqwFLje8Z4fuDDW4zkN8NZa/9nx85tKqf8CrgL+6dpby4bDwnNUFe/iSOa72G3GEHF7ay1HMt8FIDJ+hjtDE0OMzWanpKqJ+KhAlxbJPm4LSmu9A5gN/A44B9gONCul6pRS9UAzsAE4F/gVME1rvc1lkR5fGnDgmHNZGAuLXcrX28jv7R1WV99KiBMqzvmsOzl1sds6KM75zE0RiaHqs02HufO3q1i3u8il9+lzDEprbQPeA95TSgUD0zEW69qAUmCv1rqhj18xEIIwkmVPzYDLaw/5+hiTI1rbJUEJ92tvrT2p80KcqpIq4yO3pr7VpfdxdhYfjkTkiTX3moBj6wwFAI2uvrGvt5Gg2iRBCQ/g4xfWazLy8QtzQzRiKOt0zGD2srh2aexQWHibCahjzo13nHeprvVP9U3trr6VECcUP3Y5JvPRa/JMZm/ixy53U0RiqOpOUF6uTSFOt6A82BrApJS6B3gGuBxjuvn7rr5xeLAfALUNrm3mCuGMrokQxTmf0d5ai49fGPFjl8sECdHvusbdvV3cghr0CUpr3a6UWo6xDupRIA+4ZCDqBoYE+mA2m2ho7qCj0yoFY4XbRcbPkFlaKqYAACAASURBVIQkXK6r1ygkyLVl3k46QSmlwoE6xwQKj6C13gcsGuj7ms0mwoN9qaprpbq+jVjZE0oIMQzUNbYBEBbk69L7ON0+U0r9TClVAVQAyUqpl5VSzyilhnUhuq6ddIsqXD4nQwghPEJto9GCCvWEBKWUug9jV93/Btocp98FLgUed01og8PI6CAAiiVBCSGGAavN3j3uHuriLj5nW1C3AXdqrd/Asd271vpj4EbgGhfFNih0JaiicklQQoihr6KmmU6rnchQP/x8XDuNwdkElQRk93I+Hwjvv3AGn8TYYAAOl9S7ORIhhHC94oomAOKjglx+L2cT1E7g6h4/dxWhuwsY1tUoxyUaiyBzC2u7i8cKIcRQ1TXeHh8d6PJ7Ods++ymwUim1FPAFHlNKjcfYyPAcVwU3GIQG+RIT7k95TQuFZQ2MGuHyIupCCOE2eY7eooSYYJffy6kWlNZ6K5AKbMbY2sIfo7r5eK31JteFNziMSzR6OXV+jZsjEUII18ouMD7nunqPXOlklgFHAh9pra/QWl+KUSzW9W28QSBtdAQAGTmVbo5ECCFcp7W9kyOlDZhNMGZkqMvv5+w084uA3RzdnXcesFspdaYrAhtMpo6LBmBPdgV2u4xDCSGGpsNF9dhsdpLiQvDzdX0hImdbUI8BP9dad6950lqfC/wv8IQrAhtMkmKDCQ/2paahjfwyd+8+IoQQrrE3x6ggNyE5YkDu52yCGgOs6OX8CozK4cOayWRiWqrRitqRWebmaIQQwjV2HzQSVNfnnas5m6AOAhf3cn45cKT/whm85k0aAcDmjBI3RyKEEP2vpa0TfaQaswmmjBuYBOVsJ+KvgTeVUoswtn63AzMwSh3d4KLYBpUZKgYfbws6v4bK2haiwo7dQ1EIIQavjJxKOq12VFI4Qf4DU4LV2Wnm7wJnA1bgeuAqjCS1RGv9puvCGzz8fL2YOT4GgE17i90cjRBC9K+vHZ9rsyfGDtg9T2bL99XAahfGMugtnjqSzRklrNpRwEVLxrg7HCGE6BcdnVa27jOGLxZNHTlg93UqQSml/IE7gZmAN2DqeV1rfWX/hzb4zJ0UR6C/N4eK6jhcXEdKvOvXCQghhKulH6ygqbWT5BEh3QWyB4KzkySewxiH8geagaZjDgH4eFtYOt34dvHVtnw3RyOEEP1j7c5CABZNjR/Q+zrbxXcucK3W+kNXBjMUnDVnFJ9uymPVjgKuXz4B/wFYzCaEEK5S19jGpowSTCY4fVbigN7b2RZUB8ZUc3ECYxPDGD8qnKaWDlZvl1aUEGJwW7OzgE6rjZnjY4kJDxjQezuboJ4EfquUGpjJ74PcxUuNCRIfbTiETbbgEEIMUna7nc+3GEtdz547asDv72z/05XAFKBUKdUAtPe8qLWO6e/ABrP5k0YQE+5PcWUTOw6UMWdinLtDEkKIk7Yzq5zC8kYiQvyYnTZw08u7OJugnnFpFEOMxWLmwsWjeeGj/bz91UFmp8ViMplO/EYhhPAg/1mTA8DFS0bjZTmZzS/6h1MJSmv9sqsDGWrOmZfMu6uz0fk17DhQxuw0aUUJIQaPg/k1ZORWEuDnxTnzkt0Sg7ProAKAO4A0wOI4bcLYXXeG1tqlBWOVUq9gdDN29jg9RWt9SCmVBLwAzAPKgbu11p+6Mh5n+Pt68b0zxvHCR/t5dWUWsyZIK0oIMXi8s8qYF7d8fjKBA1Ta6FjOttmeBR4GYjFq74UDszFKHr3rksiONgO4RGsd1OM45Lj2JrAXY0PFH2DUDBw9ADGd0PIFKUSE+HKoqI5NUkRWCDFIHMyvYcu+Uny8LW6tiuNsgjofuE5rfRHGdPMHtdaTgBcBl06Md1SxGA+k93ItFZjliKfdUY7pI+BWV8bkLF9vC1cuSwXglU8y6ei0ujkiIYQ4sVc+zQTgwkUpRIT4uS0OZydJBAF7HM/3YySFDIzp51981yCUUj5Abztg2YHRGF17/1RKzQMKMBLSCowux3ytdc9qFlnAnO8aU385e14yK74+TGF5Ix+tP8TlZ4xzd0hCCHFcew5WsCe7kkA/Y5jCnZxtQeUBkxzPszBq8gHYgP4oOLcAKOnlKAKCgQ3AI0A8xu6+byulpmIkzuZjflczMLCryfrg7WXmB5dMBuCtrzRVdS1ujkgIIXpntdp4/qN9AFx2+jiCAnzcGo+zLahngdeVUjcBHwAblFKVwOnAzu8ahNZ6LccUoD1Gz1bae0qpm4GLgH0Y9QF7CgAav2tM/WmGimHepDi27CvlpU8y+em1M0/8JiGEGGCfbc4jr6SemIiA7oID7uTsflBPAj8EqrXWO4EfAedhtFZuc114oJS6UCl14zGnfYBWIBNIcoxTdRnvOO9Rbr1oEt5eZtbuLGSXLnd3OEIIcZTahjZeXZkFwG0XTcLX23KCd7ieUwlKKXUD8L7WegOA1volrfUc4HsYEyhcyQI8pZSao5SyKKWuxegSfEtrrTHGxh5TSvkqpU7H2Jr+dRfHdNLiIgO55mwFwDPvpNPc2uHmiIQQ4hsvrthPU0tHd4+PJzhuF59j4oIXRtfbi8A6pVTFMS+bATwO/NlVAWqtP1BK/R/wBhCHMQZ2gda6qxLr5RjbgZQDlcCtWut9rornu7jstLFsyighp6CWl1Zk8qPvTXV3SEIIwY4DZazeUeAYM5/kMWs2+xqDugHjg7+r2umh47zuk36NqBda678Cfz3OtQJguatj6A8Wi5kfXzWdHz+5ls8257FwSjxTU6X+rhDCfZpaOnjmHWMVz/XnjichJtjNEX3juF18WuvngdOAZRitqO8BZ/Q4TsdYrHu5y6McQkaNCOHqs4yuvj+9sYu6xjY3RySEGM5e+GgfVXWtqKRwLl461t3hHKXPWXxa6/UASqkUjPVG3XtHKKWitNaVLo5vSPreGePYpcvJPFzNk2/s4sFb52E2e0aTWggxfGzaW8yX2/Lx9jLzP1dPx+Jhn0POroNqBl5SSk11TFRYCZQppXKUUi6twzcUWSxm7rt+FsEB3uzMKueDdbnuDkkIMcyUVTfzl7eNrr2bzk8jMdZzuva6OJug/ooxfbsBuAZYiLEN/AbgL64JbWiLCvPnx1fPAIyyIgcOV7s5IiHEcNFptfGHV3fQ1NLBnLQ4LlzsEeVLv8XZBHUWcJujQOtlwGda6y8xqjoscFVwQ92ciXFcvGQMVpud37y8jcpaqTIhhHC9lz/JJOtIDVGhfvzP1dM9ZtbesZxNUCagWSnljTFp4jPH+SBAPlW/g5suSGPK2ChqG9p47MWttLZ3nvhNQghxitbsLOCDdblYzCbuvX4WIYHuLWfUF2cT1Hrg98A/AG/gQ6XUFOApYLWLYhsWvCxmfn7DbOIiA8gprOPpt9Kx2+0nfqMQQpyk7IIannaMO/3gkslMHB3p5oj65myCugOjosMM4CqtdTVwPcaY1N0uim3YCAn04YFb5uLva2F9ehFvfKHdHZIQYoiprm/lNy9uo6PTxjnzRnHegmR3h3RCzm75XoJRQqjnuZ+5JKJhalRcCPdeN4vHXtzKG19owkP8WD4/2d1hCSGGgObWDh55fguVda1MSI7gjkuneOy4U099lTp6AnhEa93keH5ckqz6x5yJcfzw8qn89d09PPveHkIDfVgwJd7dYQkhBrFOq43fvrydQ0V1jIgK5P9unoO3l7OdZ+7VVwtqNsZ4U9fz45EBk3507vxkahvbeG1lFn94bSePBPoweUyUu8MSQgxCdrudp99OZ/fBCkKDfHjkB/MJDfJ1d1hOO26C0lqf3ttz4XpXnZlKTX0rn27K41cvbOXR2+czPrm3DYeFEKJ3drudf364j9U7CvD1sfDgrfMYERXo7rBOSl9dfEuc/SVdJZFE/zCZTNx+6RSaWjpZt7uQB5/bLElKCOE0u93Oy59k8vGGQ3hZzPzixjmkJoW7O6yT1lcX39pjfrZjrIeyAVaM7j8b0I4HbbE+VFjMJu65Zjp2u5316UVGkrpjPuNHSZISQvTtjS80763JwWI2cf8Ns5gxPsbdIZ2SvkbKgnscNwMZwHzAV2vtC0wDtgP3uDrI4cpiMfOTa2eweNpIWto6eei5zVISSQhxXHa7nVdXHuCNLzRmE9x3/SzmThrh7rBOWV/bbTR1HcCvgB9orbdqra2O63uBu4BHBybU4cliMfPTa2ewaGo8za2dPPCPTew4UObusIQQHsZmM8ac3vryIGaziXuuncnCqYN7FrCzcw1DMBbqHisYJ9dSiVNnsZi597qZnDUnifYOK7/+11bW7ixwd1hCCA9htdr4y9u7u8ec7r9hFqfNSHB3WN+Zs8nlXeBFpdQ9wG6Msai5wJ+AV1wUm+jBYjFz95XTCAn04b01Ofzx9V3UN7Vz0ZIx7g5NCOFGHZ1W/vjaLr7eW4yPt4X/u3kOM9TgHHM6lrMJ6m7g78CHPd7TATwP3OeCuEQvTCYTN10wkZBAX15csZ9/friP8poWbr5wosdtNCaEcL36pnZ+89I29h+qIsDPi4dum0daimfX1zsZzpY6agFuUkrdDSjH6SytdaPLIhPHddnpYwkL9uHpt9P5cH0uRRWN3Hf9TAL8vE/8ZiHEkFBc0cgjz2+huLKJiBA/HrptHqNHhro7rH51UuNHWusGYIeLYhEn4YxZSUSHB/D4S9vYcaCMnz29gV/eOo/YCJnxL8RQt/9QFY+9uJWG5g5Gx4fyy1vnEhXm7+6w+t3gKMgkejV5TBR//J+lJMQEcaS0gXufWk9GbqW7wxJCuNCXW4/wwLObaGjuYNaEWB6/a+GQTE4gCWrQGxEVyO//ewnTUqOpbWzjgWc38Z812bKnlBBDTEenlWfeSecvb6fTabVxwcIUHrh5zpDu2ve4KeKOmYJLtdaX9DiXBLwAzAPKgbu11p86rpkw1mndDvgALwL3aa2Hzda0Qf7ePHzbPF5dmcW7q7N5cUUmB/Kq+fHVMwj0H7p/eYUYLsprmvnty9vJLqjF28vMjy6fwplzRrk7LJdzqgWllDqslPqVUmqCqwJRSgUppX4P/LGXy28Ce4FI4AfAm0qp0Y5rtwOXYWymOA6j8vovXBWnp7JYzNx4fhoP3DyHQD8vtuwr5Z4n15FbWOvu0IQQ38FuXc49T64ju6CWmHB/nrh78bBITuB8F99DwCxgj1Jql1Lqp0qp/l6i/AmQgrGtfDelVKrj3g9qrdu11quBj4BbHS+5Efiz1rpQa10BPIyxA/CwNHfSCJ685zRS4kMoqWri3r+s5z9rcrDZpMtPiMGko9PGix/v56F/bqa+qZ0ZKoYn7zmNsQlh7g5twDg7zfwV4BWlVDRwFXA18Bul1AbgNeA9rXV9X79DKeUD9Fbp1K61LgOu0VoXK6UeBnoWj0oD8h0ll7pkAXN6XM885lq8UirCsTX9sNM1LvWvj/bx6aY8Xlyxn51ZZfzk2hlEhg7NwVQhhpLC8gb+8NpOcgvrMJvgmnPGc+WZqcNuveNJTZLQWldorZ8BbgCeABYA/wRKlFLPORLY8SwASno5ihy/u/g47wsCmo8518w3FdSPvd71fFjPt/b1tvDDy6fyy1vnEhrkw96cSv7r92v4es/x/jMLIdzNbrfz+ZY8fvzkOnIL64iJCOC3dy3mmrPVsEtOcBKTJJRSI/mm9TQT2IZRReJNIA74K0bX2/ze3q+1XotRIulkNQHHfu0PABqPc70rMckiYmBOWhxP33s6T725m51Z5fz2le0snBrPHZdOJjzYz93hCSEcKmtb+Ou7e7qLQZ82I4E7L5syrCc6OZWglFLrMVpAecCrwLVa65weL6lSSj2DMdOuv2UCSUopf0dFC4DxfNOtl4lR3eLrHtdKtNYyO8AhPNhYZf7ppjxeWrGfr/cUsze7gh9cMpnTZiRgMg2/b2ZCeAq73c4XW4/wr4/309zaSaCfF3dcNoXTZya6OzS3c7YFlQH8XGu9uY/XrMPYI6pfaa21UmoP8JhS6n8xEuXFfNNS+zdwr1JqFUZr6mHHOdGDyWTi/IUpzJoQyzPvpJN+sII/vb6L9buL+NHlU4kOl7EpIQZaaVUTz7yTzp5sY4H93Ilx/PDyKTJW7ODsJIm7nHhNBVDxnSPq3eXAcxhroCqBW7XW+xzXngVigU0Y3XvvAA+6KI5BLzYigEdvn8+q7fk8/9F+dhwo40dPrOKqsxQXLxmDt5es3RbC1To6bXy4Ppc3v9S0tVsJCfThjksns3jaSOnR6ME03CsOKKWSgcOrVq0iIWHw759yMqrrW3nu/Qy+3mtMnBgZHcSdl01mWurQKNUvhCdKP1jOs//JoKjCGCZfPG0kd1w6mdAgXzdH5h6FhYUsW7YMIEVrndfzmsdVkhADJyLEj/tvnM1uXc4/3t9LUUUjv/zHZhZOjeeWCycSEz6sJ0IK0a8qalp44eN93TNpR0YHcvulU4bM3k2uIAlKMF3F8PS9p/PBulze+uogX+8pZtv+Ui5eMobvnTFuWM8iEuK7amnr5IO1Oby3Noe2diu+PhauOjOVS5aOwdurt43KRRdJUAIAby8LVyxLZemMBF5ekcn69CLeXZ3NF1uPcPVZinPnJ8v4lBAnwWq18cW2fF7/PIvahjYAFkwZwa0XTZLeCSdJghJHiQkP4L7vz+LipWP418f72X+oiuc+yODjjYf4/vIJLJwSj3kYLhgUwll2u51t+0t56ZNMCsuNcabUpDBuvmAik8ZEuTm6wUUSlOhValI4j/9oIVv3l/LSikyKKhp54t87GBUXzDXnjGf+pBGSqITowW63k36wgje+0BzIM6qsxUUGcMN5aSyaGi+z806BJChxXCaTiXmTRjBrQixfbcvnra8OcqS0gd++vJ2U+BCuOXs88ybFyT88MazZ7Xb2ZFfw+uffJKbgAB+uPiuV5QuSZZzpO5AEJU7Iy2Lm3PnJLJudyBdb83n7q4McLq7nNy9tY3R8KJefMZaFU+KxWGSMSgwfXYnpjS80mYe/SUyXnjaG8xemDOmNBAeKJCjhNG8vC+cvTOGsOUl8vuUI76w6yKHiOn7/6k7+HXmAS5aO5cw5Sfh6yzdGMXRZrTY27inm/XU55BbWARAc4M2lp42VxNTPJEGJk+bjbeHCxaM5Z94oVu0o4P01OZRUNfHsf/byxhdZXLhoNOfOTx62Cw/F0NTS1smXW4/w4fpcymuMsqBhQb5ctGS0JCYXkQQlTpmPt4Xl85M5e+4oNmcU897qbHIK63h1ZRZvf3WQJdMTOH9RyrDaYE0MPWXVzXy26TCfbzlCY0sHYCyyvWTpWM6YlYiP9Bi4jCQo8Z1ZzCYWTR3Jwinx7M2p5IN1uew4UMZX2/P5ans+E5IjOH9hCgumxMtaKjEo2Gx2dh8s59Ov89h+oJSuinATkiO49LSxzJ0YJ7NYB4AkKNFvTCYTU8dFM3VcNMWVjXz6dR5fbTvCgbxqDuRVE/bRPpbNSuTMOUkkxAS7O1whvqWusY01Owv4dFMeJZXGJt5eFjOLpsVz/oIUxif3tim4cBVJUMIl4qOCuO3iSVx37njW7ipkxcZD5Jc28N6aHN5bk8OE5AjOnpvEoqkj8fOVv4bCfaxWGzt1Oau257NtfymdVqO5FB3uz/L5yZw1ZxRhwTKe6g7yySBcyt/Xi+Xzkzl33iiy8mr4ctsRNqQXdbeqnvsgg0VTR7J0egKTxkYNy22thXvkl9bz1fYC1uws6C5FZDbBjPExLJ+fzOy0OPn76GaSoMSAMJlMTEiJYEJKBLddPImNe4r5cusRso7U8OW2fL7clk9EiC+LphnJalximCwAFv2utKqJDelFbEwv5lBxXff5kdFBLJudyBmzEmWzQA8iCUoMuAA/b86eO4qz546ioKyBdbsKWbe7kNKqZj5af4iP1h9iRFQgi6bGM3/yCMYmSLISp66ytoWNe4rZkF7Iwfza7vMBfl4snjaSM2cnoUaFy98xDyQJSrhVYmww1y+fwHXnjie7oJZ1uwpZn15ESWUT76zK5p1V2USF+jFv0gjmTRrBxDGReEnFCtEHu91OflkD2/aXsnV/KfpITfc1Px8LcybGsXjaSGaoGJki7uEkQQmPYDKZSE0KJzUpnFsumsS+nEo2ZRSzZV8plXWtrPj6MCu+PkyQvzez02KZnRbHtNRoggN83B268ABWq43Mw9Vs3V/Ktv2llFQ1dV/z8TIzKy2WJdMSmDkhBj8f+dgbLOT/lPA4FrOJqanRTE2N5o5Lp5BTWMuWfSVs2VdCQVkja3YWsmZnIWYTjEsMZ7qKYYaKITUpTOoBDiMllU3sPljObl3O3pxKmls7u6+FBPowJy2OORPjmJ4aLTNFByn5vyY8mtn8TcvqhvPSKCxvYOu+UnZmlXMgrwqdX4POr+HNLzWBfl5MGRfNtNRoJo6OJDEmWBZTDiF1jW3sO1RF+sEK0g+WU1rVfNT1kdFBzJ0Yx9xJcahRETIDbwiQBCUGlYSYYBLOCObyM8bR0tZJRm4lu7PK2aXLKa5sYnNGCZszSgCjgGdaSiSTxkSSlhLJmJGh0sIaRMprmtl/qKr76Nr8r0uQvzdTU6OZnhrD9NRoYiJkl9qhRhKUGLT8fb2Mbpy0OMCYQrz7YAX7cirZf7iKqrpWtjoGyo3XW1BJEYxLCmNcYhjjEsOJDPWT2VseoLWtk9yiOrILasjOr+XAkWoqHAVZu/h4Wxg/Kpwp46KYnhrDmIQwaSUNcZKgxJARFxnI8vmBLJ+fjN1up6z66G/gxZVNpGdXkJ5d0f2e8GBfxjqS1bjEMFLiQ4gIkaTlSi1tnRwpredwcT3Z+TVkF9SSX1qPzX706wL9vUlLiWBiSiQTx0QyZmSY1HIcZiRBiSHJZDIRFxlIXGQgy2YnAVBd30pWXjU5hbVk59eSXVhLTUMb2zPL2J5Z1v3eIH9vRo0IISkumOQRIYyKM57LjMGT09FppbSqmfzSBvJK6skrqSOvpP5bY0dgjDWOHhHiaN2Gk5oUxqi4EBlDHOY8LkEppe4BlmqtL+lx7gzgS6Bnm/93WutfKaVMwK+A2wEf4EXgPq11J0L0EBHix4Ip8SyYEg8Y62VKKpvILqglu6CWnMJa8krqaWzp6G519RQe7MuIqMBvjshvHoOGafKyWm1U1LZQXNlEcUUjRRWNFFc0UVzZSHl187daRWAUX02KDWbUiGDGJhgJKWVkiEz/Ft/iMX8jlFJBwEPAT4GPjrk8A3hHa311L2+9HbjM8Zo24H3gF8CjrotWDAUmk4n46CDio4NYOiMBMJJWdX0rR0oaOFJa7zgayC9toKahjZqGtu7tvXsKDvAmOjyAqFB/IkP9HMfRzwP8vAZV12Freyf1je3UNrZRVddCRU0LFbUtVDqOitoWaupbe01CYNS1i4sMICEmmJR4oyWaHB/CyOggWWwtnOIxCQr4BKgA/gGMOObaTCD9OO+7Efiz1roQQCn1MPAykqDEKTCZTI7E4s+M8THd5202OxW1LZRWNlFS1UTJMY8NzR00NNdxqKjuuL/by2IiOMCHkEAfggN9up+HBPoQ5O+Dn68FPx8Lfj5e+Pl44evj+NnXCx8vC2az0RVmMZsxm02YTV0/G0mv02rHarUZjzYbVqudTsdjW7uVlrZOWto6aW7rpNXxvOuob2qnrrGNuqZ26h2Pbe1WJ/57QWSoH/FRQcRHBzIyOoj4qEDio4OIiwyUMSPxnQxYglJK+QC9baZi11qXAddorYsdCebYBDUDiFZK/RAwAW8BD2it24A0ILPHa7OAeKVUhNb62191hTgFZrOJ2IgAYiMCmEr0Udfsdjs1DW1U1rZQVddCVV2r8by+laraVqrqWqiub6W13drdChsMvCxmQoN8CA30JTLMj+gwf6LC/IkO8zdai2H+RIT4SRISLjOQLagFwJpezlsBL611cW9vUkp5AYUYXXcvAvHAO4Ad+BkQBPQcde16HgBIghIuZzKZiAjxIyLEDwg/7uvaOqw0NLXT0NxOveOxocl43tjSQWu7lda2TlrbO2ltt9LWbjWet1lp77Ris9mx2e1YrcajzWbvPgdgsZjxMpuMR4sJs9l4tJjN+Hpb8Pfzwt+39yMk0JuQIF9CAo2EFBrkg7/v4OqSFEPPgCUorfVajNbPyb6vE1jW41SOUuox4HcYCaoJ6Fkfv2u13tGr+oRwM19vC76OVogQ4sQ8vm2ulBqplPqDo4uwiw/Q6nieCage18YDJVrrWoQQQgxanjRJ4niqgOuAZqXUo0AK8ADwL8f1fwP3KqVWYbSmHnacE0IIMYh5fAtKa90KLAeWYCSr9RhjUH9yvORZx8+bgGyMFtWDAx+pEEKI/uRxLSit9cO9nEsHTjvO620Y66cecmlgQgghBpTHt6CEEEIMT5KghBBCeCRJUEIIITySx41BuYEFoLS01N1xCCHEsNPjs9dy7DVJUI6yStddd5274xBCiOFsBJDb84QkKNgOLAZKMMouCSGEGDgWjOS0/dgLJrv9OLXyhRBCCDeSSRJCCCE8kiQoIYQQHkkSlBBCCI8kCUoIIYRHkgQlhBDCI0mCEkII4ZEkQQkhhPBIkqCEEEJ4JKkk4QJKqXuApVrrS3qcOwP4Emjp8dLfaa1/9f/tnXu0XVV1h78ooMYoLRZoUZCWgb8YSoCAgFgxUqC1IK8AFSwlEB4GlCIiCgQEQqBKyCA0IEWoPMawICCCQAO1EAxIR4nyTvghI0VATHlFE57hkf4x1zEn5+5z7iM39+4zmN8Yd9yz99qPOfc+Z8295pp7TkkjgKnAEUQ5++8DX7f95hCK3ZY2+mwEXApsDzwLfMX2LaWt1vo0kHQFsD/QLNdY2ws76VdnJG1BFPEcCywEDrXd4w39uiLpUOBfgdebVh8N/DswC9iXyPgyw/bZQy9h35G0LXCT7fXK8lp00EHS/sBZRFaFO4GJtp8dcsHbUKHPupSFBwAACzhJREFUe4ClwLKmzX5ue9fSvsr6pIEaRCSNIgonfg24saV5HHCN7S9U7HoEsE/Z5nXgeuAk4IzVJ23v9KLPVcA9wG7AXwE/lrSl7YXUVJ8KxgF72Z5d0dZJv1pSOsAbgPOICtQTgNskfdT2kmEVru+MA861/c3mlZLOBgRsAqwNzJb0G9tXDIOMHSkPaJOA6S1Np9NGB0ljiAeizwHzgG8T38GdhkzwNnTQZ3PgRdt/WrHPoOiTLr7B5Wbgz4knwFa2Bu5vs9/BwHm2n7b9HHAacORqkbB/VOoj6WPANsCptpfZvp0wYJPKJnXV5w9Ieh8wmop70gf96sp4YE3b59l+w/ZVwCPA3w+vWP2i3e/kYGCa7cW2nyA6y1p9p5o4HZgMnNmyvpMO/wD8xPZdtl8DTgQ+JWnTIZK5E+306dSnDYo+OYLqB+UJdZ2KpuW2/w84wPYzkk6jZElvYhywrqTJwAjgamCK7deBMcD8pm0fBTaQtI7tFwdbjwaroM8Y4EnbL7fIvG1T+5Dr00on/YC/IFx735O0PfAUYZBuonf96soYYEHLukeJJ93aI+ndhGvyIEkzgFeAS4gHpD+j53eqrnpdZPtUSeMbKyT9EZ11GEOMNACw/Yqkp0r7r1a7xJ3poU9hHLCepAeB9YGfAcfa/g2DpE8aqP6xA3BHxfq3gDVsP1O1k6Q1gKcJV9f3gQ2Aa4iO8gRgFPFjbND4PBJYnR36gPShp7yU5ZFt2odKn1Y66fd3wFzi6fABYA/gh5I+Se/61ZVulbvBukSndjnhIv444bJcq7S3fqdqqVeb382o8r+dDrW9dx36gZeBuwnX/RvA+UQfty2DpE8aqH5gew4x+unvfm8Cf9206nFJ0wi/7AnEjX5fU3vjJr40MEn7LNccBqAPPeWFkPmlNu1Dok8rfdDvtqbP10k6hDBUD9NZv7rS232pNbYXAZ9pWnW/pH8h5jGg53eqK/QqNEbj7XTountn+7jmZUnHAc9J2pBB0ifnoIYASR+WNL24nBqsBbxWPs8nJk8bjAZ+a/t3QyVjP5kPbFTmcRqMZoX7ovb6SPq8pINbVjfuSW/61ZXW6w7dITcAkjaTdHrL6sY9WUTP71RX6AVgezGddVjp3kkaCWxEjXWUdIakjzetavRvjd/QKuuTI6ih4QXgi8Arks4gAg+mAP9W2q8Ejpf0X8STx2llXS2xbUkPANMknUi40vYEPlk26QZ93g3MlLQA+AURSLADcJjtJ3vRr67cAYworwXMIqL4xhJul27gd8DXJD1NRIBtBRwDfJkI9vhWme8YBRwPzBwuQQfIlbTX4QfAXWWe5x7gbOA+248Nh6B9ZCywjaQDy/JM4Gbbz0kaFH1yBDUElCiWzxGhvy8Qk4nXADPKJheV5Z8TE4jzgVOHXtJ+MYGYI3iWmMieZPvh0lZ7fWz/GDiZeL9mCRFKv7vtJ8smnfSrJbaXEd+zCcRc38lEGP1zwypYHymT63sQkW1LgOuAqbavJb4/DxOG6t7SdtEwiTpQ2upg+yHg0LL8PLAZsN/wiNlnJgGLgceBJ4j3oQ6CwdMnK+omSZIktSRHUEmSJEktSQOVJEmS1JI0UEmSJEktSQOVJEmS1JI0UEmSJEktSQOVJEmS1JI0UMk7AkmjJC1vJLyUNEdSa/mAqv1GSDpc0ntXu5ADRNI2km4f4L4TJT1fPo8v12hUWV4uaffy+TJJ1w6e1NUyDMKx3iXpfyS1ZtRIupA0UMk7lX2IRLG9sSNwMTXNulIygF9MvJQ7EK4mMk/3xj8Bhw3wHEOG7beJ5KXd9hJvUkEtf3RJsrrpR9mPgSTTHUr2At62fc9Adrb9KitXeW633e8HcvzhwPZNkmZKGl8SBiddShqopCuQdBBwGbCd7XmS1ibSxlxl++sV248k0v/vR5SlntLSPgeYZ/t4SRsQNYd2JEqg/CdRZnwkK8p1LJV0iO3LSq67ycDGRHbmm4HJtl+SNJHIHfdD4KvAmsBs4MhGfSlJ+wKnAB8jyrKfZPuG0rYdkQJra6JG1feA6WVkUMVXgB816XUaK2ruTCbSz0wFfgl8l6iD9TPgi7ZfLPJOt/0nbY7fOO5lwCjb+5blXctx/5JIZTOrHGd5X65BL+f6DjCRyGz+KvC/RAqkGcBHgJ8W3aYDnwd+Cxxluzk7/fXEqG9Ob+dL6ku6+JKuwPaVwH8AF0l6F1HWfAkthqeJC4FPE3WfJhAdZTu+S9SI+gTRKW4MnEsYiAllm02AqyUdQCS/PQ7YlOhI92Tl6q5jy7l3Bg4n3IlfApC0E+FWu5IwJBcTdajGSFoPuJXozDcnEqUeTZRk6YGkDxLl6FtL1u9OJCMdR7i6ZhAG5ChgV8L4HdvhenRE0o7ALcBPiISuJxEG96imzdpeg16O/Q3ClbiL7ebii1OBA4FdiMrBDxK5HrcG7iPyJTYzG9i51GJLupS8eUk38SUi0eYVwP7A9qUi8UqUjvtAYG/bd5d1RxJZlavYmOjknrC9rBihD9h+S1LDFfis7VclPQNMLJV3AX4t6U5WnsdZEzi81Dd6RNJsoiOFePK/0XYjQGNmCUoYSRije21PLW2/KtnUzwf+uULurYiHzNYSBq8RlU3fLPWUpgCzbM8t1+IWInnnQDkGmG27UQL8MUkfIQzVBWVdp2tQiaRJ5Rg7236gpXma7XvLdnOBD9o+vyxfAOwr6QO2l5bt5xNGejQx0k66kBxBJV2D7aeBE4nSJefa/mWbTUV0kM3t84B2brJTiHIbz0u6niir8VAbGe4EnpJ0pqRrS7mO3YjyHQ2Wlo65wZIiD4Qhu7flmNNszyOMxnhJLzX+iJHBhyR9qEKc9YGXS7b8Zp4oRTJhRVXThU3trwHvqdKvj2xGT2N/F7BBKW0Ona9BFWsTI9m3iZFrK483fX6FnvrAyjq9UP6v1+GcSc1JA5V0G1sS7rjPFldfJ5oDHN4qfz2wfSOwIeEGfItwi91StW2ZX5kLrEO4HL8A3Niy2bIOsiwj5rmqWIMowbBl099YwpVYFaTwNtVBHG+02XawqAqqaMjRuCedrkE79gCeJNy3rbTq1Js+DTkq73nSHaSBSrqGMn9zCDFi2ZT280qPEh3kdk3rNqfiCb6853QO8GHbl5YggL2BXcqcUKsxORo4x/ZRti8l5kI2pe/Rfo8Rc0PNMtwq6VhgATDa9uONP6Im1alUd8iLgPe3VP4dChbQs3jjDkTtrMUDPObvbc8m5rH2l/S3qyAfwLrl/6KOWyW1Juegkq6gROVdAlxo+1ZJ3yTmb24oHfkfsL1U0iXADEmLCffShVSMXErU2RhglqRjiIi/A4gCbM8TUXoAW0v6BeE6+mzZZwRhJMfQxiVYwXnAXElfJkZguxHRg8eWcx8j6XxiLmdDIrrwhjZRfA8SI4stgP/u4/kHg3OAeZKmEAEf4wjXayOKb8AHtn23pCuACyWtyjzZFqwoppd0KTmCSrqFs4g5hkbU3iVEYMOlkqpGL18lXG8/IiLjLqfa7QRRGXQREb78IGEYditG4SHgJuA24AgidHk5Maf10yLT2bSMitpR3lf6R1aUMT+UqHq7oMyx/Q2wDfAAEel3NW0i7mwvIeZ+PtOXcw8Wtu8novL2IwIQziKi7M4apFOcAPwx8K1VOMaORCBHuvi6mKyomyRdjKT9gFNsjx1uWepCmZv8NXBgI3Ix6U5yBJUk3c11wAhJnx5uQWrEnsDCNE7dTxqoJOliihvyMMLF9o6njJ5Opg8vBSf1J118SZIkSS3JEVSSJElSS9JAJUmSJLUkDVSSJElSS9JAJUmSJLUkDVSSJElSS/4fOnptE/iw7soAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_trajectory(results.R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "R    [1993548154.0907593 meter, 53558254525.37807 m...\n",
       "V    [9671.881139414403 meter / second, 2805.422229...\n",
       "dtype: object"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error = results.last_row() - system.init"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "53595343660.133934 meter"
      ],
      "text/latex": [
       "$53595343660.133934\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "53595343660.133934 <Unit('meter')>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error.R.mag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "10070.535172486145 meter/second"
      ],
      "text/latex": [
       "$10070.535172486145\\ \\frac{\\mathrm{meter}}{\\mathrm{second}}$"
      ],
      "text/plain": [
       "10070.535172486145 <Unit('meter / second')>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error.V.mag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>The solver successfully reached the end of the...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                                                 True\n",
       "message    The solver successfully reached the end of the...\n",
       "dtype: object"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Ralston\n",
    "\n",
    "system.set(dt=system.t_end/1000)\n",
    "results, details = run_ralston(system, slope_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3ib1fXA8a/kvfdeie3kOonj7El2CIQVdlllt6WUDlpK16+lg7Z0sVraUlrKJlBGGIEEQgbZeziJ45s4juO997Yl/f54leCYDAUsS7bP53n0SH5fSe/BMTq69557r8lmsyGEEEK4G7OrAxBCCCFORxKUEEIItyQJSgghhFuSBCWEEMItebo6AFdTSvkAU4AywOLicIQQYqjxAOKAHVrrjp4nhnyCwkhOG1wdhBBCDHGzgY09D0iCMlpOvPLKK8TGxro6FiGEGFLKy8u55ZZbwP5Z3JMkKHu3XmxsLImJia6ORQghhqrPDbFIkYQQQgi3JAlKCCGEW5IEJYQQwi1JghJCCOGW3K5IQik1FViutY62/+wDNAGdPZ62WWt9kf38V4DfY9TRfwrcobWu7N+ohRBC9DW3SVBKKRNwN/CXXqfGArVa68/VgCulRgPPApcAO4E/Aq8BC5wbrRBCCGdzmwQF/Bq4DPgt8PMexycBe8/wmq8C72utNwIopX4K1CmlRmitjzgzWCGcxWq10dlloePErdO47+62YrHasNlsWG02bFaw2h8DeJhNeHiY8TSb8fAw4elh3Ht5mPH18cTX2wNfb0/MZpOL/wuFcIw7JaintdYPKaXm9To+EYhWSmUDMcB64H6tdQkwGqPlBIDWulUpVYTR6pIEJVzKZrPR1tFNXVMHdY3t1DV2UNfUTkNLJ82tnTS3ddHc1kVLaxfNbcbPbe3ddHZbnRqXr7cHvj6e+Hl74ufjSaC/F8EB3gQFeBMc4E2wv/E4JMCHsGAfIkL8CPL3wmSSxCb6l9skKK116RlOtQCbgN8AXcBfgWXAVCAQaO31/FbA30lhCnFSV7eVmoY2KmpbqaxtpbKujcq6VipqW6lpaKOuqYOOzi+2vKO3lwc+Xh74eNvvvTzw8jRjNpswm02YTGA2mTCbjMcAFqsNi9VGt8WKxWKl22LDYrXS1W2lvdNCe0e3cW+/1dNx9iB6xuNpJjzEl4gQP8KDfYkI8SU2IoDYCH9iIwKIDvPHy1NqrkTfcpsEdSZa6x/0/Fkp9QOgSimVhJG8/Hq9xB9o7qfwxCBns9mobWynpKqZkspmiiubKbY/rqprxXqODal9vD0ID/IlNMhojYQH+RIc6EOgnxeB/l7GvZ/3ycd+vp74eHk4rbVitdro6DKSVVtnN61t3TS1dtLY0vnZfYtx39DSQW1jOzUN7bS2d1Ne00p5Te/vgwazCSJC/YiLCCA2IoCkmECSYoJIigkiKtRPWl/iC3H7BKWU+g2wVGt9yH7I237fDuQAqsdz/YFk+3EhzktbRzfHyxo5VtrAsVLj/nh5E20d3ad9vtkEkaF+xIT7ExXmR0yYP9Hh/kSH+REV5k9YkA9+Pp5u9eFsNpvw8zG69sLO43VtHd3UNLRR02AkrOp6o+VYXtNCeU0L1fVtVNUZt+y86lNe6+fjQWK0kaxSYoNISwglNTGEIH/vM1xNCIPbJyggC5islLrZ/vOTwAda6yql1KvARvu41RbgEWCP1vqwa0IVA0V7Rzd5xfUcLqzjcGE9+aUNlNe0YDtNiyjI34uEqEASo4NIiA60Pw4kNiJgyHRr+fl4khgdRGJ00GnPd3Vbqao3WlilVc0UVTRRVGHc1zd3cKSoniNF9ae8Jjrcn7SEEOOWGEp6YiihQT798Z8jBoiBkKDuxhh3ysOI9wPgGwBa6/1KqbuAp4EEYBtwvYviFG7KZrNRXNlMbkEturCOw4V1HC9vwtqrf87DbCIpNojh8cEMjw9heHwww+JC5EPTAV6eZuIjA4mPDGSiij7lXENzB8WVzRRWNHGstIH8EqOFWmkfu9uy/7NFrOMiAsgYFsaoYeFkDAsnOTYYD6k6HLJMttN9ZRxClFLDgGOrV6+W1cwHCavVRmFFEweOVnPgaA0H82uobz61IMBsNjEsNpiRKWGMTAolPSmUxOigIdMicjWLxUpJVTNHSxo4WtxAXnE9R4vrae9VVOLn44lKCSMzNYKs9ChGJIfi6SH/RoNJcXExCxcuBBiutS7oeW4gtKCEOKey6hZ251aw90gVB/NraGrtOuV8WJAPo4dHoFLCGJkcRlpiCL7e8ufvKh4eZpJjg0mODWb+pCTASFoFZY3kFtRyqKCOQ8drqaxtZe/hKvYergJy8fPxYExqJFnpxm14fIjM6xrE5P9QMSC1dXSzP6+a3bqS3bmVlNW0nHI+MsSXzLRIMtMiyEyLJD4ywK2KFcTneXiYSUsMJS0xlMtmGcdqGto4VFBLdl412UeqKalqZuehCnYeqgAgyN+bSRnRTBoVw6SMaCm8GGQkQYkBo6qujW0Hy9h6oIyD+TV0Wz7rng7082L8yCgmqGiy0iOJCfeXhDQIRIT4MWtcArPGJQBGwtp3pJrsvCr2Hammur6NdbuLWbe7GLMJVEo4U0bHMHlUDMPiguVvYICTBCXcls1mjCVtPVDG1v1l5BU3nDxnfBiFMVFFMzEjmhFJYTKYPgREhPixYHISCyYnYbPZTraoduRUcDC/hkMFtRwqqOXFDw8RHe7PBVnxXJAVx8jkMElWA5AkKOF2CssbWb+nhA17Syit/qzrzsfbg4kqmumZcUweFUNwgHTnDGUmk+lk6ftVc9Npbe9iz+EqduYYXYCVta0sW5fHsnV5RIb6MTMrjguy4slICZdxqwFCEpRwCxW1razfU8z6PSUUlDWePB4c4M20MbFMz4xj3MgofLw8XBilcGf+vl72FlM8FquN3IJaNmWXsjm7lOr6Nt5bn8976/OJCPFl7oRE5k9OYlhcsKvDFmchCUq4TGt7F+v3lLBmZxGHCmpPHg/w82Lm2DjmTkgkMy0CDykrFufJw2xiTGoEY1Ij+NqSTA4X1bFpn5GsKuvaeHtdHm+vyyM1PoT5kxOZOyGRsGBfV4ctepEEJfqVzWbjYH4Nq7YXsim79ORiqt5eHkwfE8ucCQlMzIjGy1NaSqJvmM0mMlLCyUgJ564rxnCooJa1u4rZsLeE/NIG8t9r4Ln3DzJeRXPxtBSmjomVuVZuQhKU6Bf1TR2s2n6cT7YXnjKulJkWwaKpycwYG4+fj/w5CucymUyMHh7B6OERfOOqTHbkVLBmZxG7civYnWtMWQgP9uHCqSlcPC2F6HDZGMGV5BNBONXhwjqWb8xnw95Sui3GPkfhwT4snJLMhVOTiY8MdHGEYqjy8vRgZlY8M7PiaWju4NPdxazYUkBxZTP/++Qwb6w+zKSMGC6ZMYxJo2KkStQFJEGJPtfVbWHjvlI+2HgMXVgHgMkEU0fHsnhGChNVtIwrCbcSEujDkjlpXDE7lYP5NazccpxN2aUnJwXHRQRwxexULpyaLC39fiS/adFnmtu6+HDTMd7fmE99k7H2XaCfF4umpXDpzGHERgS4OEIhzs5kMtlXIInk682ZrN5RxIebj1FW08Iz7+znlY9yWTw9hctnpRIZ2nsrOtHXJEGJL62moY131+ezckvByb2ThsUFc/ms4cydmChr3okBKSTQh2vmp3Pl3DS2HSjjnU+PcqiglrfW5vHOp0eZNS6B6xeOIEVK1Z1GPjnEF1ZS1cxba46wdlfxyfGlcSMiuW7BCMaNiJKZ+2JQ8DCbTo5VHS6s491Pj7Ixu5RP9xTz6Z5iZoyN4ysXjiQ9MdTVoQ46kqDEeSurbuG1VZp1u4qw2ozxpQvGxXPt/HRGJJ3PPq1CDCwjk8N48NbJ3F5nrFLx0dbjbNlfxpb9ZUweFcMNi0aSkRLu6jAHDUlQwmGVta28tkqzemcRVqsND7OJi6Ymc+38dOKjpBpPDB3RYf7cc3UW1y8cybJ1eazYUnCyoGKiiua2S0eRJi2qL00SlDinusZ2lq7SrNp2nG6LDbMJFk5J4sZFSgofxJAWHuzL3UsyuW7BCN5df5TlG48ZW8DoSuaMT+CWSzJkKsWXIAlKnFF7ZzfvfHqUt9Ycob3TgskE8yYmcuNFigRpMQlxUkigD7ddOpqr5qbzxurDfLDpGOv3lrApu5SLp6dw4yIlSyl9AZKgxOdYrTbW7CzipRWHqG1sB2DamFhuvXQUKbFSsSTEmQQHeHP3kkyumJ3K0o80a3YW8uHmAlbvLOL6BSO4el463rLgscMkQYlTHDhazb/fOUB+qbH3UlpiCHdfkcnY9EgXRybEwBEd5s/3bpzAVfPSeOnDQ2w7WM7LK3NZtb2Qu5dkMj0zVqpcHSAJSgBQ19TOc+8fZO2uYsDYMv22y0Yzd0Ki7J0jxBeUEhvMz++aRnZeFc8s28/x8iZ+//x2xo+I4utXZZIsPRJnJQlqiLNYbazcfIyXVhyipb0bL08z1y8cyTXz02XvJSH6SFZ6FE/+YB4rtxTw8spc9h6p4juPruPKOWncfLGSyexnIL+VIexIUR3/eHPfya3UJ2VEc8/VWcRFSmWeEH3Nw8PMZbNSmTU+gVdW5vLR1gKWrctjc3Yp9103jgkq2tUhuh1JUENQV7eFpR9r3lqbh9VqIzLUj29clcn0zDjpFxfCyUICffjWdeO4cGoyf/vfXgrKGnnomS0smJzE3UsyCQ7wdnWIbsPtEpRSaiqwXGsdbf/ZG3gKuA6wAI9prR/p8fyvAL8H4oBPgTu01pX9HvgAcbiwjide20NRRRMmE1w1N42bL86QFZqF6Gcjk8N4/PtzWbYuj6Uf65P7Ut1zdRazxye4Ojy34DafSkopE3A38Jdep34NKCANCAFWKqVKtNYvKqVGA88ClwA7gT8CrwEL+i3wAaKzy8KrH+WybF0eVhskRAXwvRsmMmq4LMsihKt4ehhjvhdkxfPUG/vYf7SaP720k+055dxzdRaBfl6uDtGl3GlTnl8D9wK/7XX8duB3Wus6rXUBRgK7x37uq8D7WuuNWut24KfABUqpEf0U84BQWN7IA0+u5621eQBcPS+dJx+YL8lJCDcRHxXI7+6dybeuzcLby4N1u4r5zl/Wsj+v2tWhuZTDLSillBnIAKIxutrKgTytta2PYnlaa/2QUmpej2uGYnTd5fR4Xi4w1v54NEbLCQCtdatSqsh+/kgfxTVg2Ww2Ptp6nH+/e4DOLgvxkQF8/6aJZAyTxCSEuzGZTFwyczhZI6J49JVdHCmq5/+e3sRVc9O59ZIMvDyHXlXtOROUUmoO8D3gQiCoxykbUKeUWgn8Q2u9+csEorUuPc3hE+vptPY41gr49zjfyql6nh+ymlo7+dv/9rJlfxlgrJ13z9VZMtYkhJtLiArkT9+ZzeurDvO/1YdZti6P/Uer+cltU4gJH1ofbWfs4lNKjVBKrQH+C+QD1wAJgC9GAhgG3AGUAq8ppdYqpUb2cXwt9vueW1f6A809zvfe1rLn+SFJH6/lu4+uY8v+Mvx8PHnglkncf+NESU5CDBCeHmZuWZzBH789i+hwf/KK6rn/sXVszyl3dWj96myfWC8Dv9Faf3CG80X223Kl1I+Bq+yvmdpXwWmt65RS5RhFEiX2wxl81uWXYz8HgFLKH0jm1C7BIeWjrQU8/fZ+ui1WVHIYP/zqJFlxXIgBKiMlnCe/P5fHl+5he045Dz+7jesWjOCrizPw8HCnEgLnOFuCmu7o+JL9ecuUUu/0TVineAn4pVIqG6NL74fAk/ZzrwIb7eNWW4BHgD1a68NOiMOtdXVb+Ney/Xy09TgAl10wnLuXZOLlOfj/iIUYzAL9vfm/O6eybF0eL36Yw5trjqCP1/Hj2yYTEujj6vCc6owJ6nTJSSkVAHzuN6K1rj3Ta/rAQ8CjwEGMLslngKft19uvlLrL/nMCsA243gkxuLWahjYeeX4HurAOL08z9103joVTkl0dlhCij5jNJq5dMIKRKWH8+aWd7D9azQNPrucXd08b1DsMmGy2c+cUpdQi4N9AUu/XAzat9YAtL1FKDQOOrV69msTERFeHc97ySxr49X+2UtvYTlSYHz+7fSrpSbKTpxCDVU1DG799bjt5RfX4+Xjyo1snM3lUjKvD+sKKi4tZuHAhwHD7VKKTHO3/+QewHWMC7Mwetxn2e+ECOw9V8JO/b6C2sZ0xqRE8fv9cSU5CDHIRIX488q0LmDUunraObh5+divvrT+KI42NgcbRsq444LKhOLbjrlZsKeDpt7OxWm3MnZDI924cPyTnSQgxFPl6e/LgVyeTFKNZ+rHm3+8eoKy6ha9fNXZQbY/jaIJ6H1gESIJyMZvNxksrDvHGamMe8g0XjuSWxRmyyKsQQ4zZbOLmizNIiArkidf2sHzTMRpbOrn/pomDpjjK0QT1QyBbKXUDxpwoa8+TWuu7+jow8XlWq42nl2WzYnMBHmYT9103jkXTUlwdlhDCheZOTCQs2Iff/nc76/eW0NTayU/vmDoo5j2ezxiUF9AEBGCsKNHzJpzMYrHy5Ot7WLG5AC9PMz+7c6okJyEEYGyI+PtvXUBIoDd7Dlfx86c30dza6eqwvjRHU+wi4MIvu5yR+GK6ui38+eVdbNlfhq+3Bz+/cxrjRka5OiwhhBtJTwzlT9+ezS+e2cLhwnp+8a/NPHzPTAL9B+7+Uo62oIqBdmcGIk6vq9vKIy/sYMv+MgJ8PXn4npmSnIQQpxUfFcgfvjWLuIgA8oob+Pm/NtM0gFtSjragfgD8Ryn1W+Ao0NXzpNZ6yC4t5EwWi5VHX9nFjpwKgvy9ePiemaQlShm5EOLMosL8+P23LuBn/9zE0eIGfv7Pzfz23pkEDcCWlKMtqPeA8cCbwB7gALC/x73oY1arjb/+by+bskvx9/XkN9+Q5CSEcExkqDFXKj4ygPzSBn797620dXS7Oqzz5miCGn6aW2qPe9GHbDYbT7+dzZqdRfh6e/Crr82QCbhCiPMSEWK0pKLD/dGFdfz++e10dVtcHdZ5cTRB+Wqtj/e+YWxaeM+5XizOzysf5bJii1Gt9/O7psnOt0KILyQixI+HvzGD0EAf9h6u4tFXd2OxDpwVJxxNUJ8qpTJ7HlBKXYSxrYUkqD60attxXl91GLPZxE9un8K4EVIQIYT44uKjAvn1N2bg7+vJpn2l/OedgTMq42iCehFYp5SaqJSKVkotBVYA6zH2ZxJ9YHduJU+9uQ+Ae6/JYuroWBdHJIQYDFITQnjo7ul4ephZvukYyzfmuzokhziUoLTWPwL+DKwBDgGjgNla6zu11lVOjG/IOF7WyB9e3IHVauO6BSNYPGOYq0MSQgwiY1Ij+N4N4wH49zv72ZVb4eKIzs3hBZu01n/EKDcPAH4kk3b7TnNrJ797bjttHd3MGZ/ArZeMcnVIQohBaN6kJG5YNBKrDf744k6Olze6OqSzOuM8KKVUFXC60TQz8L5SquHEAa11tBNiGxIsVht/eWUXZTUtpCaE8J0bxg+q1YiFEO7lloszKK1qYcPeEh55fjuP3T8Xf18vV4d1WmebqPvDfotiCHtl5SF25VYS5O/Nz+6Yiq/3wF/gUQjhvkwmE9+9YTxFFU0UlDXy19f38uPbJrvljghn2/L9hf4MZCjalVvBG6uPYDbBj2+bTEy4v6tDEkIMAb7envzk9il8//FP2ZRdynsb8rlyTpqrw/qcwbFpyABU19TOE0v3AHDL4lFSTi6E6FcJUYHcf+MEAJ57/yCHC+tcHNHnSYJyAavVxhNL91Df3EFWeiTXLhjh6pCEEEPQzKx4lsxJxWK18egru2h3s+WQJEG5wPsb89mtjXGnH9w8EQ8pihBCuMjtl44mJTaI0uoW/rv8oKvDOYUkqH5WWt3Mix8eAuC7N4wnIsTPxREJIYYyby8PHrhlEp4eZlZsLmDnIfeZH+VQyZhSyg/4JjAJY2fdU77ya62/0vehDT5Wq42/vr6Xzi4L8yYlMj0zztUhCSEEw+NDuPWSDJ5bnsNTb+zlHz9a4Bal5462oJ4Bfgv4Aa1AS6+bcMCKLQUczK8hNNCHr1851tXhCCHESVfOTWdkcig1De28vDLX1eEAjm9YuBi4WWv9rjODGcxqG9t54QNjX8dvXptFcMDA2zxMCDF4eZhNfPv68dz/+Kcs35jPvImJjEwOc2lMjiaoLuCwMwM5F6XUXcC/gI4eh+8DlgJPAdcBFuAxrfUj/R/h2b3wQQ5tHd1MGR3DBVnxrg5HCCE+Z3h8CFfNSePtdXn8/Y19PPb9uS4t4nK0i+9x4A9KKVdO1pkIPKq1DuxxewH4NaCANGAKcLtS6jYXxvk5h47VsmZnEZ4eZr52Zea5XyCEEC5y00WK6DA/8ksb+GR7oUtjcTRBfQWjm69cKVWvlKrseXNifD1NAvae5vjtwO+01nVa6wLgL7jRHlVWq41n3skG4Op5acRHBro4IiGEODNfH09uv2w0YCzF5sqt4h3t4nvKqVGcg1LKA8gCblVKPYZRqPEfjC6/OIyNE0/IBdymAmHTvlLyihsID/blKwtHujocIYQ4p9njE3hvfT66sI631h7hq4tds8OCQwmq57p8SqlgwKy1rndaVJ8XBewEXgCuwdiP6l3gRKVBa4/ntgJusaidxWLl5ZXGnKebLlL4+shCsEII92cymbhryRh+/NRGlq07yqUzhxMe7NvvcTj8iamUuhf4GRBv/7kSeFJr/QcnxXaS1rocmNvj0F6l1N+AS+w/95zt6g80OzsmR6zeWURpdQtxkQFcODXZ1eEIIYTDRg+PYHpmLFsPlLNsXR53L+n/8XOHxqCUUj8E/gD8DZgNzMEonPiRUup7zgvv5PXHKKV+3euwN9AOlGMUSZyQwaldfi7RbbHy2ioNGPuveHrIoh1CiIHlxkXGR+uHmwuob+o4x7P7nqMtqPuAb2qtl/Y4tkkpdRxjAu+TfR7ZqeqBB5RSxcCzwATgu8C3gYPAL5VS2UAgxj5Wzo7nnDbuLaGqro3E6EBmj09wdThCCHHe0hJDmTo6lu055bzzaR53XD6mX6/v6Nf6KGDHaY7vAhL7LpzT01qXAEswqvMagbeAh7XWbwIPAQcwEtUO+7mnnR3T2dhsNt5amwfANfPSZYdcIcSAdeNFRnHXB5uO0dzW1a/XdrQFdQC4Hug9AfYGjKo5p9NarwEmn+Z4O0YL777+iMMRe3QVBWWNhAf7MG+S0/O3EEI4zYikMLLSI8nOq2b1jsJ+3djQ0QT1EPCBUmoGsMV+bAbG3KhrnBHYQPb+xnwArpidhpenh4ujEUKIL+fyWalk51XzwcZjXDErtd96hRzq4tNafwwsxFhm6FaMZYUagSla6+XOC2/gqaxrZVduBZ4eZhZJ5Z4QYhCYOiaW6DA/ympa2JXbf9txOFxmrrVeD6x3YiyDwsfbjmOzwQVZ8YQE+rg6HCGE+NI8zCYunTmc5z/IYeWW40wZHdsv1z1jglJK/Q/4mta60f74jGQ/KIPFamPVNmPtqotnpLg4GiGE6DsLpiTx4opD7MqtoKG5o1++gJ+ti68FsPV4fLabAA4craa2sZ24yAAyUyNcHY4QQvSZsCBfJoyMwmK1sWFvSb9c84wtKK31nad7LM7sxD/a7PEJmExSWi6EGFzmT0piV24la3cVcfmsVKdf72xdfJc6+B42rfWKPopnwLJYrGzOLgOQiblCiEFpWmYsfj4eHC6sp7K2lehw5y57erYiCUer82zAkK+l3n+0mqbWThKjA0mJDXJ1OEII0ed8vT0ZPzKaLfvL2JFTzmVObkWdrYtPFo87DzsPGdtizRgbJ917QohBa+roWLbsL2N7ToXrEpRSyuG2m9a69dzPGtx2ayNBTVDRLo5ECCGcZ/KoGEwmyM6rprW9C39fL6dd62ytpGag6Ry3E88Z0qrq2iiqaMLPx4OMlHBXhyOEEE4TGuTDyKQwui1Wco7VOvVaZxuDWsBnZebiLPYdqQIgKz0KL0/pGRVCDG6ZaRHowjoO5tcweVSM065ztjGodU676iCTe9z4FjFG5j4JIYaAMakRvLU2j4P5NU69ztnGoCqB0VrraqVUFWdpTWmth/TAiz5eB4BKCXNxJEII4XyjhkdgMsGRono6uiz4eDmnkPtsXXwP8tn40g+dcvVBoK2jm8LyRjzMJtISQ10djhBCOF2gnxdJMUEUljdxvKyRkcnO+XJ+ti6+F073WJwqv6QBqw3SEoKd9i1CCCHczfC4EArLmyhwRYLqSSkVCtwPjAE+t0Kg1npJH8c1YBRVGI3M5BiZnCuEGDpS4oJgDxSUNTrtGo5ut/EqMAn4EKh2WjQDUFGlkaCSJEEJIYaQ4fEhABSUuj5BzQEWaa23nPOZQ0xxRTMgCUoIMbQkRgcCUFbd7LRrODppJx+Q9XtOo7zG2G0kPjLAxZEIIUT/iQjxw2SC2sZ2ui1Wp1zD0RbUvcBTSqm/AceAU6Kx77Y75NhsNmob2wGIDPVzcTRCCNF/vDzNhAf7UtPQTm1Du1NWNnc0QU0ExgLPnubckF3NvK2jm/ZOC77eHvj5OPqrFEKIwSEq1I+ahnYq65yz9YajXXy/AH4PxABBvW7BfR7VAFHTYLSewoN9ZQVzIcSQExpkFHU3tHQ65f0d/drvDTyvta5yShQDVFOr8Y8SFODt4kiEEKL/BfgZK5m3tnU55f0dbUH9E7hfKeWWK6EqpcYppbYopVqUUvuVUlP647rtnRYA/Lyle08IMfQE2LfaaGnvdsr7O/rJOgK4HLhNKVUAnJIutdZT+zYsxymlvIF3gScwyuGvBT5WSqVorZ1XoA90dBr/KD7eQ3IITriRmtLdlOatoLO9Hm/fUOLTLyEifqKrwxKD3Im9oFqc1IJyNEHtt9/c0TzAS2v9hP3n15RS3wZuAP7tzAufaEH5SgtKuFBN6W6O57yJzWp8SHS213M8500ASVLCqby9jE61rm6LU97foU9WrfWvnXL1vjEaONTrWC5G1aFT2ezru5vcsuNTDBWleStOJqcTbNYuSvNWSIISTnViJ/GPtx3njsvH9Pn7n/GjVSn1iVLK4YGQUX0AACAASURBVK47pdQFSqk1fRPWeQkEem853wr0fc1jL2azUblntcq+jsJ1Otvrz+u4EH2lvNpYqKCptf+7+H4C/EspZQHeBlYCOVrrTgCllA8wDpgLfNX+mq87JcqzawF6z5L1x9iO3qk87AnKIglKuJC3b+hpk5G3r2z/IpxrRlY872/IZ9HUZKe8/xlbUFrrncAU4I/AxcAOoFUp1aCUasRopWwAFgMPA+O11tudEuXZ5QCq17EM+3GnOpmgnLTMhxCOiE+/BJPZ65RjJrMX8emXuCgiMVSc2GIozklLvZ11DEprbQXeAt5SSgUBEzAm61qBciBba910lrfoD2sBk1Lq+8BTGFV8WcAyZ1/Y39f49bU6qcRSCEecGGeSKj7R3zrtxRFens4ZiHe4/MyeiNxuzT2tdadS6hLgaeA3QAFwVX9MKg7yNybonpiwK4SrRMRPlIQk+l2TfQWJQD/nLFYwKOqjtdYHgFn9fd2TCcpJy3wIIYQ7q2/qAD5b8qivSYH0l3BiiaOGlk5sNimUEEIMLXUnElSgJCi34+fjSZC/F13d1pPfJIQQYqhwuxaUUirMXdfkc4UY+xLzFbW9p2IJIcTg1dbRTX1zB54eZsKCfZ1yDYcTjVLqR0qpKqAKGKaUekEp9ZRSyutcrx3MYsKN8spySVBCiCGktMqYahoXGXByyk1fcyhBKaUexNhV97vAib6sN4GrgUecEtkAcaL+v6TS6fOChRDCbZTYE1RidKDTruFoC+prwDe11kuxb/eutX4fuB24yUmxDQjD4oz9Go+VNrg4EiGE6D/F9i/l8U6apAuOJ6hk4MhpjhcCYX0XzsCTmhACSIISQgwtecXG8lonPgOdwdEEtQu4scfPJ2qq7wN292lEA0x8VCDeXh5U1rXJhF0hxJBgs9k4UmgkqJHJzmujOJqgHgAeVEp9BPgAv1NK7QbuAn7krOAGAg+zidR4o5tPH69zcTRCCOF8VXVt1Dd3EOTvfbKS2RkcSlBa623ASGALxu61fhirm2dorTc7LboBIjMtEoADR6tdHIkQQjhf7vFaAEYkh2IyOaeCD85vHlQE8J7W+nqt9dUYi8U6b3RsAMlMiwDgwNEaF0cihBDOt/ewsdTpWPuXc2dxtMx8CbAHY9uNEy4F9iilLnRGYAPJqGHhmM0mjhTX09runI27hBDCHdhsNvbYE9SEkVFOvZajLajfAT/WWp+c86S1Xgz8FPiTMwIbSPx9vVDJYVitn/3DCSHEYFRS1Ux1fRshgd4Mj3deBR84nqDSgOWnOb4cY3PAIW/amFgAth0oc3EkQgjhPDsPVQIwLj0Ks5NWkDjB0QR1GLjyNMcvAY73XTgD17RMI0HtPFQhO+wKIQatTftKAJg+Ns7p13J0P6jfAq8ppWZhbP1uAyZiLHV0m5NiG1ASo4NIiAqgpKqFA0drGOfkvlkhhOhvVXVt5B6vw9vLgymjYpx+PUfLzN8ELgIswFeBGzCS1Byt9WvOC29gmTU+AYC1u4tcHIkQQvS9TdmlAEwZFYOvj/P3uz2fLd/XAGucGMuAt2BSEq+vOszm7FK+eXVWv/wDCiFEf1m7y/jyPWt8fL9cz6FPUKWUH/BNYBLgBZwyMqa1/krfhzbwxEcFkpESRu7xOrYeKGPepCRXhySEEH0ir6ie/JIGgvy9mDo6tl+u6WiRxDMY41B+QCvQ0usm7BZMNpLSyq1SOyKEGDxWbi0AYMHkZLy9PPrlmo72QS0GbtZav+vMYAaDuRMTeW55DgfzazhW2uD0eQJCCOFsbR3drN9TDMDF01P67bqOtqC6MErNxTn4+3px4dRkAN7fkO/iaIQQ4stbte04bR0WRg8PJykmqN+u62iCehz4g1JKaqcdcNkFwwH4dHcxDc0d53i2EEK4r26LlXfWHwXg6nnp/XptR7v4vgJkAeVKqSbglI2PtNbRfR3YQJYQFcjkUTHsPFTBu+uPctulo10dkhBCfCEb95VSVddGQlRgvxVHnOBognrKqVEMQjdcOJKdhypYvvEYV89LJ8jf29UhCSHEebFabby91thM/Zr56U5f2qg3hxKU1voFZwdyNkqpFzFacd09DmdprfOVUsnAs8B0oBL4jtb6QxeEeYqMYeGMHxnF3sNVvLc+n1sWy5KFQoiBZVN2KcdKGwkP9mX+pMR+v76j86D8gXuA0cCJ+kITxu66E7XWzv70nQhcpbVeeZpzr2FspHgZMAt4Ryk1Xmvt8gqFGxcpI0FtOMoVs1MJDpBWlBBiYLBYrLyy8hAAN12k8PLsn9Lynhwtknga+BUQg7H2XhgwBWPJozedEpmdfZJwBrD3NOdGApOBh7TWnfbVLt4D7nZmTI4akxrBhJFRtLZ38/oq7epwhBDCYZ/sKKKkqoW4yICTlcn9zdExqMuAW7TWy5VSORgJYb9S6j/Al14uQSnlDYSf5pQNSMXo2vu3Umo6UGS//nKMFl2h1rrnZOFcYOqXjamv3HnFGPY+to4PNh3jsguGEx8V6OqQhBDirNo6uln6cS4At1ycgafH+Wy+3nccvWogsM/++CBGqwWM8vO+2FF3JlB2mlsJEARsAH4NxGNsnvg/pdQ4e1ytvd6rFfDvg5j6xPD4EC6ckozFauP5D3JcHY4QQpzTG6sPU9PQTnpiCLPti2C7gqMtqAIgE6P1kouxJt9zgBX40kslaK3X0Wt9v14+7vH4LaXUncAS4ADG8ks9+QPNXzamvnTL4gzW7y1hy/4ydudWMjFDqvKFEO6ppKqZZevyALjnmqx+r9zr6XzGoF5VSl0JvAPcpZT6lf34LifFBoBS6gql1O29DnsD7UAOkGwfpzohw37cbUSE+HHzRQqAf769j44ui4sjEkKIz7PZbPz7nf10W2xcOCWZjJTTjbz0H0f3g3ocuBeo1VrvAr4FXIrRnfY154UHGFWDTyqlpiqlPJRSN2N0Cb6utdYYXY+/U0r5KKXmY+z8+6qTYzpvS+akMSwumPKaVimYEEK4pfV7StiVW0mArye3XTbK1eE4lqCUUrcBy7TWGwC01s9rracC12EUUDiN1vod4P+ApUAj8ABwuda60P6Ua4FRGHOg/gPcrbU+4MyYvghPDzP3XTcOkwneXpvHsdIGV4ckhBAn1TW1869l+wGjuCssyNfFEZ1lDMpeWeeJMTb0HPCpUqqq19MmAo8ATzgtQkBr/Xfg72c4VwRc4szr95WMYeFcMmMYH24u4LFXd/PY/XNcMrdACCF6e/rtbJpaOxk/IoqLpvXfiuVnc7YW1G0YxQaN9p/zgaZet3XAKifGN+jccfkY4iICKChr5JWVua4ORwghWL+nmM3ZZfj5ePCdr4zHZHJdYURPZ0xQWuv/APOAhRitqOuABT1u8zEm617r9CgHET8fT35w80TMJnh7XR4H82tcHZIQYggrr2nh728as4juvCKT6HC3maVz9jJzrfV6AKXUcIwJsbYT55RSkVrraifHNyhlDAvn2gUjeGP1Ef7y8k6e+ME8QgJ9XB2WEGKI6bZY+fPLO2lt72bG2DgW9+NmhI5wtMy8FXheKTXOXkm3EqhQSuUppWQV1C/g5oszyEgJo7qhncde3Y3Vajv3i4QQog+9vOIQhwvriQz1c6uuvRMcTVB/x5hf1ATcBFyAsQ38BuCvzgltcPP0MPOjW6cQ5O/Nbl3JG2tkw2IhRP/Zsr+Mt9bmYTabePCrk9xySyBHE9Qi4Gv2FcKvAVZorVdhLDs001nBDXZRYX48cMtETCZ4dWUuu3MrXR2SEGIIOF7eyONLjTUWbr1kFKOHR7g4otNzNEGZgFallBdG0cQK+/FAoM0ZgQ0VkzJiuOFChdUGf3ppB0UVTa4OSQgxiDW3dvK757bT1mFhzvgErp3fv9u4nw9HE9R64M/AvwAv4F2lVBbwJLDGSbENGTddpJiZFUdLezcP/3cbTa2drg5JCDEIdVus/OmlnZRVt5CaEMJ3bnC/caeeHE1Q92AsOTQRuEFrXQt8FWNM6jtOim3IMJtNfP/GiaQmhFBW3cIfXthBV7es1yeE6Ds2m41/vLmPPYerCAn05v/umIqvt6PrhbuGo1u+l2Gscdfz2I+cEtEQ5evjyS/umsYPnviU7LxqHn11Nw9+dTIeLlxJWAgxeCz9WLNqeyHeXh784q5pbjXf6UzOttTRn4Bfa61b7I/PSJJV34gM9eNXX5/BT/+xkU37SgkOyObea7LcugkuhHB/H209ztKPNWYT/Pi2ySgXr1LuqLO1oKZgjDedeHwmMoGnD6UmhPDzu6bxy2e2sGJzASEBPtyyWKaaCSG+mA17S/jHm3sBuPfacUwdHeviiBx3xgSltZ5/usfC+camRfKjWyfzyPPbeW2VJsDPi6vmprk6LCHEALP1QBmPvrILq81YHGDxjGGuDum8nK2Lb46jb3JiSSTRd6ZnxvHt68fz1//t5dn3jN1DJEkJIRy181AFf3xxBxarjWvnp3PjopGuDum8na2Lb12vn20Y86GsgAWj+88KdGJssy762KJpKXRbjcqbZ987gNVq4xo3nrMghHAPu3IreOT57XRbbFwxO5XbLxs9IMeyz1ZmHtTjdiewH5gB+GitfYDxwA7g+84Ocii7ZMYwvn39OACeW36Qt9YccXFEQgh3tjm7lN/+dxud3VYWzxjG16/MHJDJCc4+BtVy4rFS6mHgWq31jh7ns5VS92GsKvEvp0Y5xF08fRhmk4m/vbGX5z/Ioa2zm1suzhiwf3RCCOdYu6uIJ17bg9VqY8mcVL62ZOAmJ3BwHhQQjDFRt7eg83gP8SUsmpaC2Wzir6/v4fVVh6lv6uDea7Lw8HB0rrUQYjD7YNMx/rUsG5sNblg0clB8iXU0ubwJPKeU+j6wB2MsahrwGPCik2ITvSyckkygnxd/emknH209Tn1TBw/eOhkfL9k2Xoihymq18eKHOby1Ng+AOy4bzbULRrg4qr7h6Nfv7wDbgHeBUqAEeB2je+9B54QmTmdaZhwPf3MmgX5ebDtYzi+e3ixr9wkxRHV2WfjLK7t4a20eHmYT37thwqBJTuBggtJat2mt7wAiMVpO04BIrfW3tdbdToxPnMbo4RH88duziAzx5VBBLQ88uV5WQRdiiGls6eShZ7awYW8Jfj6ePPS16Vw4NdnVYfWp8xo/0lo3ATudFIs4D8mxwfz5u3N4+L/byC9p4IEn1/OjWyczeVSMq0MTQjjZsdIGfvvcdiprWwkP9uVXX5/O8PgQV4fV52SEfQCLDPXjj/fN4oJx8bR1dPObZ7fy9toj2Gyy+pQQg9WGPSX88K8bqKxtJT0plEe/N2dQJieQBDXg+fp48uNbJ3PL4gxsNnhueQ5/eWUXbR3S8yrEYGKx2nh++UH+9PJOOrssLJicxB/vm0VkqJ+rQ3MatysRt1cKztVaX9XjWDLwLDAdqAS+o7X+0H7OBDwMfAPwBp4DHhxKY2Mmk4kbFymSY4J4fOlu1u8p4WhxAz+5fQrD4oJdHZ4Q4kuqa2rnsVd3s/dwFWazibuXjOGKWakDvoz8XBxqQSmljimlHlZKjXJWIEqpQKXUn4FHT3P6NSAbiAC+DrymlEq1n/sGcA3GZoojMFZe/5mz4nRnM7Pieez+uSTHBlFS1cwDT67nk+3HXR2WEOJL2KMr+e6j69h7uIrgAG8evmcGS2anDfrkBI538f0SmAzsU0rtVko9oJSK7+NYPgCG02tVCqXUSPu1H9Jad2qt1wDvAXfbn3I78ITWulhrXQX8CmMH4CEpKSaIR783hwunJNPZZeHJ1/fy+NLdtLZ3uTo0IcR56LZYeeGDHH757y3UN3UwNi2Svz4wj6z0KFeH1m8c3VH3ReBFpVQUcANwI/B7pdQG4BXgLa1149neQynlDZxulyyb1roCuElrXaqU+hUQ1+P8aKCw59JLQC4wtcf5nF7n4pVS4fat6YccX29PvnfjBMakRvDPt7NZs7OIg/k1fP+miYxJjXB1eEKIcyivaeHRV3aRe7wOswluXpzB9QtHDrkdts+rSEJrXaW1fgq4DfgTMBP4N1CmlHrGnsDOZCZQdppbif29S8/wukCgtdexVj5bQb33+ROPh/wK6xdOTebx++eQGh9CRW0rP/3HRp5ffpCubourQxNCnIbVauPDzcf4zl/Wknu8jsgQX37/rVncuEgNueQE51EkoZRK4LPW0yRgO8YqEq8BscDfMbreZpzu9VrrdRhLJJ2vFqB3mYo/0HyG8ycSUzOC5Nhg/vK9OSz9OJe31hzhrbV57Mqt5Ac3Txy0palCDESVta389X972HekGoDZ4xP45jVZBAd4uzgy13EoQSml1mO0gAqAl4GbtdZ5PZ5So5R6CqPSrq/lAMlKKT+tdZv9WAafdevlAArY1ONcmda63gmxDEhenmZuu3Q0U0bF8vjS3RSUNfKDJz7l2vkj+MqFI/GWtfyEcBmbzcbH2wp59r0DtHV0ExzgzbeuHccF4/p6mH/gcbQFtR/4sdZ6y1me8ynGHlF9SmutlVL7gN8ppX6KkSiv5LOW2kvAD5VSqzFaU7+yHxO9jBoezpMPzOP55Qf5cHMBr39ymA17S/j29eMZmx7p6vCEGHKKKpr451vZ7D9qtJpmZsVx7zXjCA3ycXFk7sHRIon7HHhOFVD1pSM6vWuBZzDmQFUDd2utD9jPPQ3EAJsxuvfeAB5yUhwDnp+PJ/deO465ExN56o29FFU087N/bmLR1GTuvGIMQf5DtztBiP7S3tnN/z45zLJ1eXRbbIQEevONq8Yye3zCkCgfd5RpqC+Lo5QaBhxbvXo1iYmJrg6nX3V1W3hrbR6vrzpMt8VKcIA3t14yikXTUobkgKwQ/WFHTjlPL9tPZa1Rz7V4xjBuu3TUkP1yWFxczMKFCwGGa60Lep5zu5UkRP/x8vTgxkWKWePi+fub+zhwtIa/v7mPFZsL+PpVmWSmSbefEH2luLKJ597PYXtOOQDD44P51nXjyEg53ewbAZKgBJAYHcTv772AjftK+e/7B8kvbeCn/9jErHHx3HnFGKLDhnzFvhBfWENzB699rFmxpQCL1Yafjwc3XzyKK2YNlx2xz0ESlACM9fxmj09gyugYlq07yptrjrBxXynbD5Zz+axUrls4Ysh2QQjxRXR1W1i+8Rivr9K0tHdjNsHF01O4ZXEGYUG+rg5vQJAEJU7h6+3JTRcpFk5J4oXlOazfW8Lb6/L4aGsB18wfwZLZqfj6yJ+NEGdisVhZu6uY11ZpKuzjTONHRnH3kkxZvPk8ySeNOK3oMH8evHUyV81L48UPD7H3cBUvrTjE+xvzuXGR4qJpKXh5SveEECdYrTY27C1h6ce5lFQZK7MlxQRx1xVjmJQRLdV5X4AkKHFWI5LCePiemew7XMULH+ZwpKiep9/O5s01R7hufjqLpqXIRF8xpNlsNrbsL+OVj3IpLG8CIC4igJsuVsyZkCgVsV+CJCjhkHEjo3h0xJxT/kd8etl+Xv/kMNfMT2fx9GHS9SeGlG6LlQ17S3hrzRGO2xNTVJgfNy5SLJichKcUQHxp8okiHGYymZiZFc/0zDi2Hijj9U8Ok1/SwLPvHeSN1UdYMieVS2cOl2IKMai1d3azalshyz7No6rOWH0tPNiXrywcwUXTU/DylB6FviIJSpw3s9lIVDPGxrErt5LXVmn08TpeXpHL/z45wsIpSVw5J42EqEBXhypEn2lo7uDDzQUs35hPY0snAAlRgVw7P515kxIlMTmBJCjxhZlMJiaPimFSRjTZR6p5+9M8dudWsmJzASs2FzBldAxXzkkjKz1SBojFgHWkqI7lG4+xfk8J3RYrACOSQrl+4QimjYnDLGNMTiMJSnxpJpOJcSOjGDcyisLyRt7bkM/anUXsyKlgR04FKbFBLJ4xjPmTkgjw83J1uEKcU1e3lU3ZpSzfmI8+XgeAyYR86epnshbfEF6Lz5kamjtYsaWADzYdo76pAwBvLw/mTkhg8YxhjEgKlf/Bhdsprmzik+2FrN5ZdPLvNsDXk0XTUrh05nDiIgNcHOHgI2vxiX4XEujDjYsU184fwdYDZazcUkB2XjWrtheyanshqfEhXDQtmdkTEof0hmzC9do7utm4r5RV24+Tc6z25PHk2CAun5XK/ImJUqHqIvJbF07l5Wlm9vgEZo9PoKSqmZVbCli9o4j80gaeXraf/7x3gEkZMcyfnMTU0TEy0Cz6hdVq4+CxGj7dXcz6PSW0dXQD4OvtwaxxCSyalsyoYeHSyncxSVCi3yREBXL3kkxuvWQUW/aXsWZXEXt1JdsOlrPtYDmBfl7MGp/AvImJjBoWLoPPok/ZbDbyiutZv6eEDXtLqGloP3kuIyWMRdNSmDUuHn9fGSd1F5KgRL/z9vJg7sRE5k5MpLaxnfV7ilm7s5j80gZWbilg5ZYCwoN9mDE2nguy4hmdGiGz8cUXYrPZOF7exKZ9pazfU0xpdcvJc9Fhfswen8CCyUkkx8oaee5IEpRwqfBgX66am85Vc9MpKGtk7c4iNmaXUlnbygebjvHBpmOEBvowY2wcM8bGkZkWKWsAirOyWG3kFtSy9UAZWw+UUV7TevJcaJAPs8bFM3dCIiolTLrw3JwkKOE2hsUFc+cVY7jj8tEcLW5gU3Ypm7JLKatuYcWWAlZsKcDPx4PxI6OZMiqGyaNiCAuWbQsEtLZ3kZ1XzXZ7d/GJibQAIYHeTB0dy+zxCWSlR8oeTAOIJCjhdkwmE+lJoaQnhXLbpaMoKGs8uTdVQVkjW/aXsWV/GQDpiSFMHhXL+JFRqJQwWf9siLBabRwtqWe3rmSPriK3oBaL9bMpM3ERAUwfG8f0zFhUSrh0EQ9QkqCEWzOZTAyPD2F4fAi3XjKKyrpWdh4yJgBnH6kir7iBvOIGXlul8fX2YExqBFnpUYwbEcnw+BAptBgkbDYbZdUt7D9aQ3ZeFXsPV53SSjKbYNSwcCZlRDN9bBzJMUHSfTcISIISA0p0mD+XzhzOpTOH097Zzf68anblVrLvSBXFlc3syq1kV24lAEH+XmSmRTJqWDijhoeTlhAq41cDhNVqo6iyiQNHaziYX8PB/GpqGztOeU5UmB8TVTQTVDTjRkQRKKuUDDqSoMSA5evtyZTRsUwZHQtATUMb2XnV7DtSxb4j1VTXt53SHejlaWZEUqiRsIaFMzI5TMaw3ERDcwdHiuo5XFhnv9XT1Np5ynNCAr0ZkxrBmNQIJoyMJjE6UFpJg5wkKDFoRIT4MX9SEvMnJZ3sEso5VkPOsVpyj9dSVNFMzrHaU1YLCA/2JS0xhPTEUNISQkhPCiU82Fc++JyovqmDgrIGjpU2kldsJKWelXYnhAf7kpkWQWZaJJmpEZKQhiBJUGJQMplMxEcFEh8VyIVTUwBobOlEH6/lUIFxO1rcQG1jO7U57ezIqTj52tBAH5Jjg0iOCSI5NoikmCCSY4NlSabz1NreRWlVC4UVTRSUNVJQ2sCxssaTa9z15OPtQXpiKCOSQlEpYYxICiM6zE8S0hAnCUoMGcEB3qd0CVqtNsprWjha3MDRknryius5WtxAfXMH9XkdZOdVn/L60EAfkmKCiI3wJy4ygNiIAOIiAoiNDBiy4x/tnd1U1bVRUdtKcWUzJVXNlNjvaxvbT/saPx9PhsUFMyw+mLSEEEYmh5EcEyTl3+Jz3C5BKaW+D8zVWl/V49gCYBXQ1uOpf9RaP6yUMgEPA98AvIHngAe11t39GLYYgMzmz1pZsyckAEa1WFVdG4UVTRSWN1FY0UhRRRNFFU1G4mruYP/Rz79XkL8XMeH+RIT4ERHiS2SocR8R4kdkqB9hQT74+XgOmBaBzWajtb2b+uYO6hrbqW/uoLahncq6NirrWqmqa6Wyru2USrrevDzNxEUGkBgdyPD4ECMpxQUTE+4/YH4PwrXcJkEppQKBXwIPAO/1Oj0ReENrfeNpXvoN4Br7czqAZcDPgN84L1oxWJlMJqLD/YkO92fyqJiTx61WG9X1bRRXNVNe00JZtXErr2mhrKaVptYumlqNkvcz8fQwExzgTXCANyGB3gQH+BAc4E2gnxd+Pp74+Xri5+OJr7cn/vaffbw98PQw22+mzx57mjGbTNhsNqw2G9jAarNhf4jFYqWjy0JXt5XOLgsdXRY6uyx0dlnp6LTQ3NZJS1sXzW1dp963dlHX1E59Uwed3dZz/r48PcxEhfkRE+ZPfFQACdGBJEQZt6gwf5l/JL4Ut0lQwAdAFfAvIK7XuUnA3jO87nbgCa11MYBS6lfAC0iCEn3IbP4scfVms9moa+qgsq6Vmvp2ahraqG4w7msa2qmub6O+uYOOTosx5nWGri934+PtQViQD2FBvoQG+RAW5ENUmD/RYX5Ehxm/i9BAH5lrJpym3xKUUsobCD/NKZvWugK4SWtdak8wvRPURCBKKXUvYAJeB36ute4ARgM5PZ6bC8QrpcK11rUI4WQmk4nwYF/Cg30h5czPa+/spqmli8aWDhpbOmlo6aSxuYOWti7aOi20dXTT1t5t3Hd009bRRUeXhW6LjW6LFYvFStf/t3f+0XaNZx7/3PzSxG1NtXQapTpd+k2jQsOgZkSq6BiKNpiSGiF+NLQGVVNEpNIwRmSJSTVVSmWtGYoqRUN/CPFj1iRt/WriISujpGT8iiYRRHLv/PG8R/Y9d59z7z335ty9V57PWmeds9/33fs8z9n7vM9+n/fZ77O+nQ1tbaxf38aGtnYGDGihBWgZ0ILbCX8fMKCFIYMHMmTwQLYYPJAhgwd02G4dOpgthw5+733LzHbFIA2NHEhBP9PMK3Af4P6c8g3AIDN7MW8nSYOA5bjr7npgOHAL7sk4F2gFsjGqlc/DgDBQQWF43xB3323zwaH9LUoQlIKmGSgzm4+Pfnq633rgC5mipZKmA5fhBupNIPuPr/hg1jQmaRAEQVAE//C96gAADPdJREFUCh/XKWk7STOSi7DCEKDiyF8MKFM3AnjJzN5oloxBEARB31MGJ/NrwHhgraSLgU8Ak4Efp/q5wDmSfoOPpqamsiAIgqDEFH4EZWZvAwcDY3Bj9SA+BzUzNZmTth8BnsVHVFOaL2kQBEHQlxRuBGVmU3PKHgPG1mjfhj8/ddEmFSwIgiBoKoUfQQVBEASbJ2GggiAIgkJSOBdfPzAQYMWKFf0tRxAEwWZHpu8dWF0XBiqtWjF+/Pj+liMIgmBz5qNAh6WYw0DBQmBf4CV8VYsgCIKgeQzEjdPC6oqW9vb25osTBEEQBF0QQRJBEARBIQkDFQRBEBSSMFBBEARBIQkDFQRBEBSSMFBBEARBIQkDFQRBEBSSMFBBEARBIQkDFQRBEBSSWEmiF0g6C9jPzI7IlO0P/Ap4K9P0MjObJqkFmAacgmcFvh74dkprXwhq6LQDcB2wN/Ay8E0zuyfVFV6nLJJuBI4GsvKNMrNl9fQsOpJ2xXOjjQKWASeaWacn84uOpBOBHwLvZIpPB/4LmA0cia/4MtPMLm2+hN1H0p7AXWa2bdoeQh0dJB0NXIKvqvAAMMHMXm664F2Qo9cWwGpgXabZI2Z2UKpvWK8wUA0gqRXPP/Ut4M6q6tHALWb21ZxdTwG+ktq8A9wOnA9cvOmk7R5d6HQT8ChwCPD3wM8l7WZmyyiwTjUYDRxhZvNy6urpWVhSx3cHcCWe2HMccJ+kj5vZqn4VrueMBq4ws+9kCyVdCgj4JLAVME/Sn83sxn6QsS7ppm0iMKOq6rvU0EHSSPzm6GBgEXAZfj3u3zTBu6COXrsAr5vZX+fs0yu9wsXXGHfjqed/mFO3O/BYjf2OB640s+Vm9gqenv7UTSJhz8nVSdKngD2AKWa2zsx+ixuwialJkXXqgKShwAhyzk839CwyY4HBZnalmb1rZjcBfwT+qX/Faoha/5/jgelmttLMnsM7yUJeZ7ghmgR8r6q8ng5fA35hZg+lLOLnAX8naacmydwdaulVr8/rlV4xgsoh3ZFunVPVbmb/BxxjZi9KmkpaDT3DaGAbSZOAFuBmYLKZvQOMxFPSV3gaGC5pazN7va/1yNILnUYCz5vZm5myp4E9M/X9olMe9fQE/gZ37f1I0t7AC7hBuouu9SwyI4ElVWVP43e2pUHSQNxFeZykmcBa4Fr8pumjdL7OiqrfHDObImlspUDSX1Ffh5H4CAMAM1sr6YVU/+wml7h7dNIrMRrYVtITwEeAB4EzzezP9FKvMFD57APcn1O+ARhkZi/m7SRpELAcd3NdDwwHbsE7x3OBVvxPV6HyeRiwqTvzhnSis8yk7WE16pupUx719PxHYAF+J/g4cBjwU0mfo2s9i0yZZc+yDd6Z/QR3G38ad10OSfXV11kh9avxX2pN77V0KPw5rNNHvAk8jLv13wWuwvvAPemlXmGgcjCz+fjop6f7rQe+kClaKmk67nc9Fz+RQzP1lZO0pjFJeyTbfBrQic4yg8u9pkZ903TKoxt63pf5fJukE3BD9RT19SwyXZ2jUmBmK4D9MkWPSfoPfP4COl9nZdKvMjKvpUNpz6GZnZ3dlnQ28Iqk7emlXjEH1YdI2k7SjORmqjAEeDt9XoxPklYYAbxkZm80S8YGWAzskOZvKoxgo6uiNDpJ+pKk46uKK+enKz2LTPU5gPLI/h6Sdpb03ariyvlZQefrrDT6mdlK6uvQ4RxKGgbsQAl0lHSxpE9niir9X+V/1bBeMYLqW14DxgNrJV2MBx1MBn6c6ucC50j6DX5nMTWVFRYzM0mPA9MlnYe70A4HPpealEmngcAsSUuA3+FBBPsAJ5nZ813oWWTuB1rSIwKz8Si+UbibpUy8AXxL0nI88uuzwBnAN/Cgj4vSPEcrcA4wq78EbZC51NbhP4GH0vzOo8ClwB/M7Jn+ELSHjAL2kHRs2p4F3G1mr0jqlV4xgupDUpTKwXio72v4ZOEtwMzUZE7afgSfIFwMTGm+pD1mHD4f8DI+aT3RzJ5KdaXRycx+DlyAP1OzCg+pP9TMnk9N6ulZWMxsHX7djcPn/S7AQ+lf6VfBekiaVD8Mj2xbBdwGTDOzW/Fr6incUC1MdXP6SdRGqamDmT0JnJi2XwV2Bo7qHzF7zERgJbAUeA5/Huo46L1ekVE3CIIgKCQxggqCIAgKSRioIAiCoJCEgQqCIAgKSRioIAiCoJCEgQqCIAgKSRioIAiCoJCEgQo2CyS1SmqvLHQpab6k6rQBefu1SDpZ0vs2uZANImkPSb9tcN8Jkl5Nn8em36g1bbdLOjR9vkHSrX0ndb4MfXCsAZL+R1L1yhpBCQkDFWyufAVfNLYrxgDXUNBVV9IK4NfgD+c2ws34itNd8S/ASQ1+R9MwszZ80dKyPcQb5FDIP10QbGp6kAqkkQV2m8kRQJuZPdrIzmb2Fh2zP9dq95dGjt8fmNldkmZJGpsWDw5KShiooBRIOg64AdjLzBZJ2gpfNuYmM/t2Tvth+LL/R+HpqCdX1c8HFpnZOZKG4zmHxuCpUX6FpxkfxsbUHaslnWBmN6Q17yYBO+KrMt8NTDKzNZIm4GvH/RQ4CxgMzANOreSaknQkcCHwKTw9+/lmdkeq2wtfGmt3PF/Vj4AZaWSQxzeBn2X0msrGXDuT8GVnpgG/B36A58R6EBhvZq8neWeY2YdrHL9y3BuAVjM7Mm0flI77GXwJm9npOO3d+Q26+K5/BybgK5u/BfwvvgTSTOBjwK+TbjOALwEvAaeZWXal+tvxUd/8rr4vKC7h4gtKgZnNBX4JzJE0AE9vvooqw5PhamBfPAfUOLyjrMUP8HxRf4t3ijsCV+AGYlxq80ngZknH4Aving3shHekh9Mxu+uo9N0HACfj7sSvA0jaH3erzcUNyTV4TqqRkrYF7sU7813whVJPx1O1dELSB/DU9NXp6w/FFyMdjbu6ZuIG5DTgINz4nVnn96iLpDHAPcAv8AVdz8cN7mmZZjV/gy6O/a+4K/FAM8smYZwGHAsciGcQfgJf/3F34A/42olZ5gEHpBxtQUmJkxeUia/jC23eCBwN7J0yFXcgddzHAl82s4dT2an4asp57Ih3cs+Z2bpkhN5vZhskVVyBL5vZW5JeBCakLLwAf5L0AB3ncQYDJ6f8Rn+UNA/vSMHv/O80s0qAxqwUlDAMN0YLzWxaqns2rax+FfBvOXJ/Fr/JrE5d8Dae0XR9yqc0GZhtZgvSb3EPvmhno5wBzDOzSurvZyR9DDdU309l9X6DXCRNTMc4wMwer6qebmYLU7sFwAfM7Kq0/X3gSEnvN7PVqf1i3EiPwEfaQQmJEVRQGsxsOXAentLkCjP7fY2mwjvIbP0ioJab7EI89carkm7HU2w8WUOGB4AXJH1P0q0pdccheCqPCqtTx1xhVZIH3JAtrDrmdDNbhBuNsZLWVF74yOBDkj6UI85HgDfTKvpZnkvJM2FjNtNlmfq3gS3y9OsmO9PZ2D8EDE+pzaH+b5DHVvhItg0fuVazNPN5LZ31gY46vZbet63znUHBCQMVlI3dcHfc55Orrx7ZAIcN6dUJM7sT2B53A27A3WL35LVN8ysLgK1xl+NXgTurmq2rI8s6fJ4rj0F4CobdMq9RuCsxL0ihjfwgjndrtO0r8oIqKnJUzkm936AWhwHP4+7baqp16kqfihy55zwoB2GggtKQ5m9OwEcsO1F7XulpvIPcK1O2Czl38Ok5p8uB7czsuhQE8GXgwDQnVG1MTgcuN7PTzOw6fC5kJ7of7fcMPjeUleFeSWcCS4ARZra08sLzU00hv0NeAWxZlQW4GSyhcyLHffA8WisbPOZfzGwePo91tKR/6IV8ANuk9xV1WwWFJuagglKQovKuBa42s3slfQefv7kjdeTvYWarJV0LzJS0EncvXU3OyCVFnY0EZks6A4/4OwZPvPYqHqUHsLuk3+Guo8+nfVpwIzmSGi7BHK4EFkj6Bj4COwSPHjwzffcZkq7C53K2x6ML76gRxfcEPrLYFfjvbn5/X3A5sEjSZDzgYzTueq1E8TV8YDN7WNKNwNWSejNPtisbk+gFJSVGUEFZuASfY6hE7V2LBzZcJylv9HIW7nr7GR4Z9xPy3U7gGUFX4OHLT+CG4ZBkFJ4E7gLuA07BQ5fb8TmtXyeZLqVqVFSL9LzSP7MxjfmJePbbJWmO7YvAHsDjeKTfzdSIuDOzVfjcz37d+e6+wswew6PyjsIDEC7Bo+wu6aOvOBf4IHBRL44xBg/kCBdfiYmMukFQYiQdBVxoZqP6W5aikOYm/wQcW4lcDMpJjKCCoNzcBrRI2re/BSkQhwPLwjiVnzBQQVBikhvyJNzFttmTRk8X0I2HgoPiEy6+IAiCoJDECCoIgiAoJGGggiAIgkISBioIgiAoJGGggiAIgkISBioIgiAoJP8PWX4NEE3SANwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_trajectory(results.R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "R    [-3723637.7203674316 meter, -1060434466.552573...\n",
       "V    [-214.86614651603597 meter / second, 0.7785432...\n",
       "dtype: object"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error = results.last_row() - system.init"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "1060441004.1725632 meter"
      ],
      "text/latex": [
       "$1060441004.1725632\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "1060441004.1725632 <Unit('meter')>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error.R.mag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "214.8675569931363 meter/second"
      ],
      "text/latex": [
       "$214.8675569931363\\ \\frac{\\mathrm{meter}}{\\mathrm{second}}$"
      ],
      "text/plain": [
       "214.8675569931363 <Unit('meter / second')>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error.V.mag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = results.R.extract('x') * 1.1\n",
    "ys = results.R.extract('y') * 1.1\n",
    "\n",
    "def draw_func(state, t):\n",
    "    set_xlim(xs)\n",
    "    set_ylim(ys)\n",
    "    x, y = state.R\n",
    "    plot(x, y, 'b.')\n",
    "    plot(0, 0, 'yo')\n",
    "    decorate(xlabel='x distance (million km)',\n",
    "             ylabel='y distance (million km)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZgcVdn+8e9ACBADIksERAQR7hgFNAEURDZBjYAgO0R4EVA2QSKLoCg7sr3IEjEgCKJgFBDCIpsXhkX5vRCWsAQeNlk0hFUgJMRAMr8/Tg10Ot09lZmu6Zr0/bmuudJ96nT1U1OZfvqcOnVOR2dnJ2ZmZmWzUKsDMDMzq8UJyszMSskJyszMSskJyszMSskJyszMSskJyszMSqntE5Sk9SS93IPXjZZ0TZ1tn5L0pqTBvY/QzKw9DWh1AK0iqQPYGzhjPl83GDgGOBS4tsb27YFzgSWbEKaZWdtq2wQFHAdsCZwIHN1VKGkx4GRgR2ARYDzww4iYnlW5AXgFOB9YoXKHkr4HHAYcD/yq4PjNzBZo7dzFNzYiRgATq8pPA9bJfgQMAc6p2L5rROwAvFRjn9cCQ4G/Nj9cM7P20rYJKiKmVJdl3X77AIdHxEsR8SZwJLCnpEXrva5in1MjYk5RMZuZtZN27uKrZTlgceAWSZWTFL4LfAJ4oiVRmZm1ISeoub0GzALWi4gAyFpOnwSebmVgZmbtpm27+GqJiNnA74FTJS2TJaczgOtaG5mZWftxgprXIcC/gEmkgRBrACOz5GVmZn2kw+tBmZlZGbXdNais225d4EXArSIzs9ZamHRP6b0R8d/KDW2XoEjJ6c5WB2FmZnP5MnBXZUE7JqgXAS677DKWX375VsdiZtbWpk6dyqhRoyD7bK7UjglqNsDyyy/PSiut1OpYzMwsmeeSi0fxmZlZKTlBmZlZKTlBmZlZKZXuGpSk9YDrI2JIne2bAbcC71QUnxoRJ/RFfGZm1jdKk6DmYwHB4cAVEbFL8VGZmVmrlKmL7zhgf9ICgo2MAB4sPhwzM2ulMiWoegsIVhsObC7pOUnPSzq9a60mMzNbcJQmQTVaCLCLpAGkiVyvBj4NbAZsDvj6k5nZAqY016DyiIj3gK9UFD0l6STgVOCI1kRlZmZFKE0LKg9JH5N0hqSBFcUDgZmtisnMzIrRr1pQpBVvRwEzJB0PrAocDfympVGZmVnTlb4FJWmUpLcBImImMBLYiJSs7gCuAM5sXYRmZlaE0rWgImICsFTF88uAyyqePwhs0ueBmZlZnyp9C8rMzNqTE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZXSgLwVJS0EDAWGALOBqcBTEdFZUGxmZtbGuk1QkjYCfgBsDixRsakT+I+km4DzIuIfzQhI0nrA9RExpM72lYGLgC8CLwMHRcRfmvHeZmZWHnW7+CStLuk24DfAM8B2wMeAxYBBwCrAnsAUYJykv0lao6eBSOqQtA9wCzCwQdVxwEPAMsB3s/f+ZE/f18zMyqlRC+r3wPERcUOd7S9kP9dL+hGwbfaa9XoYy3HAlsCJwNG1KmQJcB1gi4iYBdwm6Vpgb+AnPXxfMzMroUaDJL7YIDnNJSI6I+Jq4Au9iGVsRIwAJjaoMwx4PiKmV5Q9DqzZi/c1M7MSqtuCqjX4QdKHgEVr1H293mvyiogpOaoNBmZUlc0gdTmamdkCJNcoPklbAL8GPl61qYM0WGLhJsdVz3Rg8aqyQcDbffT+ZmbWR/IOMz8PuAf4H2BmceF0azKwsqTFI+KdrGxoVm5mZguQvAlqBWDLiHiiyGC6ExEhaRJwkqSjgA2AbYD1WxmXmZk1X96ZJK4DtigykHokjZJU2YW3PfBp0j1QFwJ7R8QjrYjNzMyKk7cFdRjwkKSdSfdEzancGBF7NSugiJgALFXx/DLgsornLwAjm/V+ZmZWTvNzDWoRYBrwoeLCMTMzS/ImqC2AzZs1nZGZmVl38l6D+hetHb1nZmZtJm8L6ofAhZJOBJ4G3q3cGBEe5m1mZk2VN0Fdm/17ZUVZJ31/o66ZmbWJvAlq1UKjMDMzq5I3QS0WEVFdKGlR4Bjgx02NyszM2l7eQRK3S/psZYGkr5KmGNq36VGZmVnby9uCuhSYkCWlfwFnAztl5UcUFJuZmbWxXAkqIo6Q9BpwGzCbtFDhl31flJmZFSVvFx8RcSppuPmHgCOcnMzMrEh1W1CSXiENIa+2EHCdpDe7CiJiSAGxmZlZG2vUxXdYn0VhZmZWpdGS77/ty0DMzMwq5b4GZWZm1pecoMzMrJScoMzMrJScoMzMrJRy3agraXFgP2AEaWXdjsrtEbFT80MzM7N2lneqowuA7YCbgLeKC8fMzCzJm6C+DuwWEeOLDMbMzKxL3mtQ7wJPFBmImZlZpbwJ6hfAKZKWKzIYMzOzLnm7+HYC1gKmSpoGzKrc2Iy5+CStDYzN3ucZYK+IuLdGvc2AW4F3KopPjYgTehuDmZmVR94ENabIICQNBMYDZwEbAdsDt0j6RERUD8oYDlwREbsUGZOZmbVW3vWg3p+XT9KSwEIR8UYT49gEWCQizsqej5P0fWBn4NdVdUcADzbxvc3MrITytqCQtD/wY2DF7PnLwNkRcUoT4hgGPFZV9jiwZo26w4Hlsng6gD8CR0fEf5sQh5mZlUSuQRKSDgNOAc4FvkzqhvsFcISkHzQhjsHAjKqyGcCgqjgGkJacvxr4NLAZsDng609mZguYvC2oA4H9IuIPFWV/l/QccCJwdi/jmA4sXlU2CHi7siAi3gO+UlH0lKSTgFOBI3oZg5mZlUjeYebLAfOMqAPuA1ZqQhyTAVWVDc3K3yfpY5LOyAZVdBkIzGxCDGZmViJ5W1CPADsCP68q35l0rai3/gZ0SBpNGjG4PWm4+dVV9V4DRgEzJB0PrAocDfymCTGYmVmJ5E1QPwNukLQ+cHdWtj5pCqTtehtERMySNJJ0H9TxwLPAthHxiqRRwPkRMTgiZmb1ziIlqxmkeQLP7G0MZmZWLnmHmd8i6SvAQcDupJtkHwPWjYhJzQgkIh4BNqxRfhlwWcXzB0nD0s3MbAGWe5h5RNwB3FFgLGZmZu+rm6Ak/QnYJyLeyh7X5fWgzMys2Rq1oKYDnRWPzczM+kzdBBUR36n12MzMrC806uL7Rs59dEbEjU2Kx8zMDGjcxXd9zn10Ags3IRYzM7P3NeriyzvLhJmZWdM16uIbVG9btYionujVzMysVxp18b3NB6P46unAXXxmZlaARglqM7pPUGZmZoVodA1qQh/GYWZmNpdG16BeBoZFxKuSXqFBayoihhQRnJmZta9GXXyHA9Oyx4f1QSxmZmbva9TF99taj83MzPpCrtnMJS0FHAJ8Bli0entEfLPJcZmZWZvLu9zG5cAI4C/Aq8WFY2ZmluRNUBsBW0TE3d3WNDMza4K80xk9Q7op18zMrE/kbUHtD4yRdC7wT2BO5cZstV0zM7OmyZughgNrAhfV2OapjszMrOnydvH9FDgZ+CiwRNXPksWEZmZm7SxvC2ogcElEvFJkMGZmZl3yJqhfAYdIGh0Rc7qt3QOS1gbGAmuRBmXsFRH31qi3Mqmr8YvAy8BBEfGXImIyM7PGHn/2dR5++lXWXG1Zhq6ydFP3nTdBrQ5sBewh6Vng3cqNEbFeb4KQNBAYD5xFGtK+PXCLpE9ExFtV1ccBdwNbAhsC10j6XEQ805sYzPqr16bcz5SnbmTWzDcYuNhSrPipkSyz4vBWh2Vt4PFnX+cnY//Oe+/NYcCAhThpvy81NUnlTVAPZz9F2QRYJCLOyp6Pk/R9YGfg112VJK0BrEO6J2sWcJuka4G9gZ8UGJ9ZKb025X6em3wlnXPSd8ZZM9/guclXAjhJWeEefvpV3ntvDnM64b335vDw06/2fYKKiOOa9o61DQMeqyp7nDRysLre8xExvaper1pwZv3VlKdufD85demc8y5TnrrRCcoKt+ZqyzJgwELvt6DWXG3Zpu6/0XIbfwV+HBH35NmRpC8BJ0TEZj2IYzBQvWz8DKB62fm89czawqyZb8xXuVkzDV1laU7a70stuQZ1JHC+pNnAn4GbgMlZ1xqSFgXWBjYGvp295rs9jGM6sHhV2SDSsvM9qWfWFgYutlTNZDRwsaVaEI21o6GrLN30xNSl7n1QETERWBc4FfgacC8wQ9Kbkt4itVzuBL4OnAB8Lm9rq4bJgKrKhmbl1fVWlrR4N/XM2sKKnxpJx0KLzFXWsdAirPipkS2KyKx5Gl6DyoaUXwVcJWkJ4POkm3XnAFOBhyJiWoNd5PU3oEPSaGAMaRTfWsDVVfGEpEnASZKOAjYAtgHWb0IMZv1O13Umj+KzBVHeUXxkiaiQOfciYpakkaT7oI4HngW2jYhXJI0Czo+IwVn17YELSPdAvQrsHRGPFBGXWX+wzIrDnZBsgZQ7QRUtSzIb1ii/DLis4vkLgPsvzMwWcHnn4jMzM+tTTlBmZlZK852gJH1EkhObmZkVKvc1KElHAIcDHwHWkHQMMA0YHRHvNnyxmZnZfMrVEpJ0OGlV3YOB/2bFVwLfAn5eTGhmZtbO8nbV7QPsFxF/IFvuPSKuA/4H2LWg2MzMrI3lTVArA0/WKH+e1OVnZmbWVHkT1H3ALhXPO7N/DwTub2pEZmZm5B8kcShwk6SNgUVJUw0NJS1k+LWigjMzs/aVqwUVEf8HrEFayXY8aUbxm4ChEfGP4sIzM7N2NT/3My0DXBsRO0bEt0iTxX6omLDMzKzd5R1m/k3gAebuzvsG8ICkzYsIzMzM2lveFtRJwI8i4v17niLi68BRwGlFBGZmZu0tb4JaDbi+Rvn1pAUDzczMmipvgnqCtDBgtZHAc80Lx8zMLMk7zPxEYJykDUlLv3cCw0lTHe1RUGxmZtbG8g4zvxL4KjAb+DawMylJbRQR44oLz8zM2tX8LPl+G3BbgbGYmZm9L1eCkrQ4sB8wAlgE6KjcHhE7NT80MzNrZ3lbUBcA25Fmj3iruHDMzMySvAnq68BuETG+yGDMzMy65B1m/i5pqLmZmVmfyJugfgGcImm5IoMxMzPrkreLbydgLWCqpGnArMqNETGkt4FI2gk4GVgBuB3YMyJerlP3eODIqji2iogJvY3DzMzKIW+CGlNkEJKGAReRZqaYCJwKjAM2q/OS4cDBETG2yLjMzKx1ciWoiPhtwXF8G7guIu4CkHQU8B9Jq0dEraXmR5BmtzAzswVU3vugBgH7AsOAhbPiDtLqusMjotsJYyUNBJausakz2+/EroKImCHpBWBNYK4EJWkFYHngSEnrA68Bp0fExXmOxczM+oe8gyTGAscCHyXNvfcRYF3SlEdX5tzHBsCLNX7+DQwGZlTVnwEMqrGfIaRrVGOAlYD9gbMkbZkzDjMz6wfyXoPaEhgVEddLmgz8LCIelnQh8PE8O8gGMHTU2iapaxn5SoOAt2vsZxKwSUXR7ZJ+R7qR+IY8sZiZWfnlbUENBiZljx8F1ske/wJoxoq6kwF1Pcm6FFfOyuciaUNJh1QVDwRmNiEOMzMribwtqGeBzwIvAI+TBilcDMwBPtyEOC4H7pK0CXA38HPggYiodXPwO6R7sp4EbiSN9NuN+iP+zMysH5qfa1CXS9oGuAbYS9KxWfl9vQ0iIh4G9sr29yrwGWDHru2Sxkq6Mat7H7A7aan5acC5pHum7ultHGZmVh55h5n/QtKLwOsRcZ+kA4ADSCPoDm5GIBFxFXBVnW37VT2/AriiGe9rZmbllKsFJWkP4OqIuBMgIi6JiPWAHUgDKMzMzJqqbgsqu29pAGnk3cWk0XKvVFUbTrpedFZhEZqZWVtq1MW3B2kdqM7s+TN16nlot5mZNV3dBBURF0p6gtQNeBuwPfB6RZVO0n1KDxcaoZmZtaWGgyQi4g4ASasCz0dEV2sKSctGxKsFx2dmZm0q731QM4BLJJ0JPELq1ttC0j9Jy1w8XlSAZmbWnvLeB/VLYCjpvqNdgS+RloG/EzinmNDMzKyd5U1QWwD7RMQzpDnvboyIW4GTSJPAmpmZNVXeBNUBzJC0CPAV0hRDkOboe6eIwMzMrL3lvQZ1B3A68AawCDBe0lrA2aQRfmZmZk2VtwW1L2mhwuHAzhHxOmkV3GnAQQXFZmZmbSzvXHwvAttUlR1RSERmZmY0nuroNOC4iJiePa7LycrMzJqtUQtqXdL1pq7H9XQ22GZmZtYjjaY62rTWYzMzs77QqItvo7w76ZoSyczMrFkadfFNqHreSbofag4wm9T9NweYBQwqIjgzM2tfjYaZL1Hx8x3SrOXrA4tGxKLA54B7gdFFB2lmZu2n0TWo6V2PJZ0AbB8R91Zsf0jSgaRZJc4vNEozM2s7eW/UXZJ0o261Jcg/G4WZmVlueZPLlcDFkkYDD5CuRX0BOBO4tKDYzMysjeVNUAcBvwLGV7zmXeBC4PAC4jIzszaXd6qjd4A9JR0EKCt+PCLebnZAWStt44jYtkGdTUkT1a4GTAJ2j4inmx2LmZm1znxdP4qIacDEIgKRNBg4BjgUuLZBvWWBa4C9snqHADdLWiMi5hQRm5mZ9b28gyT6wg3AqnQ/InA74NGIuCoi3o2I04FFSetUmZnZAqLPRuBJGggsXWNTZ0S8BOwaEVMkHQus0GBXw4DJVWUBrAnc2oxYzcys9fpyiPgGwN9qlM8GBkTElJz7GQzMqCqbgWezMDNboORKUJL+CfweuDwiHuvJG0XEBNLw9N6azrzJaBDQ9AEbZmbWOnmvQR0DrANMknS/pEMlrVhgXI1M5oORhF2GMm+3n5mZ9WN5h5lfClwqaTlgZ2AX4GRJdwKXAVdFxFvFhTmXq4HTJO2UPf4BadLaCX30/mZm1gfmaxRfRLwSEWOAPYDTSNeVfg28KOmCLIE1naSxkm7MYngZ2Bo4Cngd2AHYOiJmFfHeZmbWGrkHSUj6GB+0nkYA95BmkRgHLA/8knRf0vq9CSgijq1Rtl/V8zuAz/fmfczMrNzyDpK4g9RaepY0WGK3iHiqosprksYAFzU9QjMza0t5W1APAz+KiLsb1LmdtEaUmZlZr+UdJHFgjjqvAK/0OiIzMzPKNdWRmZnZ+5ygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslJygzMyslAa0OoBqkkYDG0fEtg3qHA8cCcyqKN4qIiYUHJ6ZmfWR0iQoSYOBY4BDgWu7qT4cODgixhYemJmZtUSZuvhuAFYFzs9RdwTwYLHhmJlZK/VZC0rSQGDpGps6I+IlYNeImCLpWGCFBvtZAVgeOFLS+sBrwOkRcXEBYZuZWYv0ZQtqA+DFGj//BoiIKTn3MwS4HRgDrATsD5wlactmB2xmZq3TZy2obABDRxP2MwnYpKLodkm/A7YjdROamdkCoEzXoHKRtKGkQ6qKBwIzWxGPmZkVozSj+ObDO8Apkp4EbgQ2A3bL/jUzswVEv2hBSRor6UaAiLgP2B04DZgGnAvsGRH3tDBEMzNrstK1oCLi2Bpl+1U9vwK4oq9iMjOzvtcvWlBmZtZ+nKDMzKyUStfF1wcWBpg6dWqr4zAza3sVn8ULV29rxwS1AsCoUaNaHYeZmX1gBeDpyoJ2TFD3Al8mzWIxu8WxmJm1u4VJyene6g0dnZ2dfR+OmZlZNzxIwszMSskJyszMSskJyszMSskJyszMSskJyszMSskJyszMSskJyszMSskJyszMSqkdZ5JoCkmjgY0jYtsGdY4HjgRmVRRvFRETCg6vx3Ie16bA2cBqwCRg94h4ul79VpK0E3Ay6U7120lrh71cp26pz5ektYGxwFrAM8BeETHP3feSVgYuAr4IvAwcFBF/6ctY58d8HNdmwK2kRUu7nBoRJ/RJoD0kaT3g+ogYUmd7vzpfXXIcV6/PlxPUfJI0GDgGOBS4tpvqw4GDI2Js4YH1Ut7jkrQscA2wV1bvEOBmSWtExJy+iDUvScNIf/gjgYnAqcA46q++XNrzJWkgMB44C9gI2B64RdInIuKtqurjgLuBLYENgWskfS4inunLmPOYz+MaDlwREbv0cZg9IqkD2Bs4o5uq/eZ8wXwdV6/Pl7v45t8NwKrA+TnqjgAeLDacpsl7XNsBj0bEVRHxbkScDiwKfKXoAHvg28B1EXFXRMwEjgK+JGn1OvXLfL42ARaJiLOy3/s44FFg58pKktYA1gF+FhGzIuI20heJvfs64Jw2IcdxZcp8fmo5DtgfOLFehX54viDHcWV6fb7cgqqSfaNbusamzoh4Cdg1IqZIOpZsZvQ6+1kBWB44UtL6wGvA6RFxcQFhd6tZxwUMAyZXlQWwJqk536caHRcp1oldBRExQ9ILpFifrNpPqc5XDcOAx6rKHicdS3W95yNielW99QqMrTfyHhekb+TLSdof6AD+CBwdEf8tNsQeGxsRP5O0SYM6/e18Qb7jgiacL7eg5rUBaabz6p9/A0TElJz7GUK65jEGWIn0jeMsSVs2O+CcmnVcg4EZVWUzgEHNCXO+NTqu+Ym1bOerWt5jKdv56U6ueCUNAP4FXA18mtRNuzlQ2utPOf+m+tv5ynVczTpfbkFVyS6IdzRhP5NI3Rddbpf0O1IX2Q293X8P4plAE44LmM68fzyDgLebsO/51ui4JI0HFq8qrhlr2c5XDdPJdyx565VFrngj4j3m7kZ+StJJpOuKRxQaYbH62/nKpVnnyy2ogkjaUNIhVcUDgZmtiKeJJgOqKhvKvN1+ZTBXrJIGAStTI9Z+cL7y/t4nAytLWrybemWR67gkfUzSGVmXbpcynZ+e6m/nK5dmnS+3oIrzDnCKpCeBG0lN3N2oP4Ksv7gaOC0bvn018ANgDjChlUHVcTlwV9ZXfjfwc+CBiHiiRt2yn6+/AR3ZbQBjSKPd1iKdg/dFREiaBJwk6ShSF+g2wPp9HG9euY6LdE1wFDAjux1gVeBo4Dd9GGvT9cPzlVdTzpdbUE0kaaykGwEi4j5gd+A0YBpwLukenHtaGGKPVB3Xy8DWpBFxrwM7AFtHxKwGu2iJiHiYNBx+LPAq8Blgx67t/el8Zb/fkaQP8NeBnwDbRsQrkkZJquwS2p7U7/8ycCGwd0Q80tcx55H3uLJRmCNJQ9FfA+4ArgDObEngvdCfz1cjRZwvr6hrZmal5BaUmZmVkhOUmZmVkhOUmZmVkhOUmZmVkhOUmZmVkhOUmZnVJGk9STWXp+nmdaMlXVNn26ckvZmtoNCQE5Qt0CQNltTZNbGlpAmSulsmAEkdkr4rabHCg+whSetIuq2Hr91T0qvZ402y39Hg7HmnpK2yx5dIurJ5UdeOoQn7WkjSPZKqZ6WwHsj+/+8D3EKaASLv6wZLOh343zrbtyfdE7Vknv15JglrN9sB7+aotxFwAfCHYsPpGUkLk+I7sIe7+COQZ1G8H9CcORwLFRFzshkLxgKbtjqeBcBxpPWpTiTNAAFA9oXtZNIN74uQ1vL6YcVs7DcAr5CW7ZlrVQRJ3wMOA44HfpUnCLegrK1ExOsRMS1H1bJ/KG8LzImIu3vy4oh4p97KwlX13oyIN3ryHn0tIq4nzWu3SatjWQCMjYgRVCxXkzmNtH7VOqQ5FIcA51Rs3zUidgBeqrHPa0nzDP41bxBuQVmpSdoduAT4QkRMlPRh4BFgXEQcXqP+INIfzI6kKYuOrto+AZgYEYdJWpH0TW8j0vpRt5JaJINIc8QBTJP0nYi4JJsvbn9gFdJs0zcA+0fE25L2BL4P/AkYTfp2eROwb9e3S0k7AD8F1iAtbf7jiBifbfsCaRqYEcALwK+BMxqsUnwQ8OeK4zqWD9a52p+0bP0JwP2kb6ufJHWtjIqI17N4z4iIZevsv2u/lwCDsw8dJH012+9nSdNHjcn205nnd9DNe50G7AlsTJob8Z/AN7Pfy0qkD7b9SSu5bk1aVuWAiLilYjdd80NO6O79rL5aS2pkK+nuA2yarSGHpCOByZIOiIj/NlqKIyKmZq/JHYdbUFZqEfE70uStYyUtRFoa/C2qEk+F84AvA98gzXE2usHufwXMBtYlfSiuQuo7fyF7LcBqwB8l7QocC/wQWJ30QboNsG/F/tbK3ntz4Luk7sT9ACRtRupW+x0pkVwA/EnSMElDgJtJH+ZrAgeTEmXNZQkkLUlaGvymqk1bkdYXGk7q6jqTlEAOAL5KSn7VM7bnJmkjUrfgdcDngR+TEu4BFdXq/g662fePSB9+W0RE5QKGJ5Am7d2CtBzKQ8A/smN5gDR3XaWbgM2z9YisuZYjLQ1yi6Q3JL0B/B+py/wTRbyhT6L1B/uRlgG/FNgJ+GKtVTmzD+7dgG9FxN+zsn1JM5nXsgrpQ+7ZiJiVJaElImK2pNezOi9HxDuSppAmj70+K39O0u2kFVG7LAJ8N/um+Kikm0gfpJC++V8bEV0DNM7OBiUMIiWjeyOiazG3J7OZrc8BTqkR9+dJXy6rl2SYCRwSEe9JOpeUxMdExJ3Z7+IvpAlze+pg4KaI6Frq+wlJK5ES1S+zska/g5ok7Z3tY/NsXa5KJ0XEvVm9O4ElI+Kc7PkvgR0kLVHRbTuZlKSHklra1jyvkVrm60VEAEhalNQ6f7qIN3QLykovIv5Fmj19FPC/EXF/naoifUBWbp9IWg6klp8COwOvSrqatMTBw3ViuB14QdKJkq6U9BjpIvLCFdWmdXVjZN7K4oGUyO6t2udJETGRlDQ2kfR21w+pZbCMpGVqhPNRYHo2Y3SlZ7OF4uCDVVqfqdg+E1i01vHl9BnmTfZ3AStKWip73uh3UMuHSS3ZOaSWa7WnKh7PYN7jgbmP6bXs3yEN3tN6ICJmA78HTpW0TJacziC1qAvhBGX9xckMh2IAAALWSURBVOdI3XGbZl19jVQOcJid/cwjIq4FPk7qBpxN6harObItu75yJ7A0qctxF9JF30q1lhzpqNhWb+mAAcBVpGPs+lmL1JX4Zo36c6g9iKPW6MR6ybkn3qlR1hVH1zlp9Duo55vA86Tu22rVx9Td8XTFUfOcW68dQlrKfRJpIMQawMgseTWdu/is9LLrN98htVguJyWUWvdZPE76gPwCHyx4tyY1vsFnF3xPA34fERcBF2UDAG7OrglVJ5MDgdMj4mcVr1+deUc51fME6dpQZQw3k5LdY6Q/8qcqtm1NGuixZ419TQU+JGnxiKiVNIryGPMupLcBaR2j//Rwn29GxE2SpgF3SrokIqqvrc2P5bJ/pzasZblExARgqYrn00gDYb7fzeuObbDtKXKOknWCslLLRuVdCJwXETdno4bOljS+8gMd0h+PpAuBMyX9h9S9dB41Wi7ZqLNhwBhJB5NG/O0KPEsanda1oNwISfeRuo42zV7TQUqSw6jTJVjDWaQP4O+TktKWpNGDh2TvfbCkc0jXcj5OGl04vs4ovodILYu1gf+X8/2b4XRgoqSjSQM+hpO6XrtG8fV4xxHxd0mXAudJ6s11srVJyfKp7ipa+bmLz8ruZNI1hq5RexeSBjZclLViqo0mdb39mTQy7rfU7nYC2Jv0TfuvpA/9jwNbZknhYeB60p303yMNXe4ktZj+msX0c6paRfVk9yvtQfrm+Shppd9tI+Kx7Brb10j3lkwijfT7I3VG3EXEW6RrPxvnee9miYgHSaPydiQNQDiZNMru5Ca9xRHAR4BjerGPjUgDOdzFtwDwirpm/ZCkHYGfRsRarY6lLLJrk88Bu3WNXLT+zS0os/7pKqBD0pdbHUiJbAM84+S04HCCMuuHsm7IfUhdbG0vaz39hBw3BVv/4S4+MzMrJbegzMyslJygzMyslJygzMyslJygzMyslJygzMyslP4/61z1ix7HJVoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "system.set(dt=system.t_end/100)\n",
    "results, details = run_ralston(system, slope_func)\n",
    "animate(results, draw_func)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
